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The Correction of Flexible Plate Supersonic 
Nozzle Contours by Influence Methods’ 


WILLIAM N. MacDERMOTT?t 
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SUMMAR\ 


Methods of obtaining the contour corrections necessary for the 
production of a satisfactory flow uniformity at the exit of a flex 
ible plate supersonic nozzle are reviewed Basic requirements 
of the methods are considered, and in particular it is found that 
contour corrections, small in comparison with the boundary 
layer thickness, produce flow changes which cannot be ade- 
quately described by the simple theory of characteristics A 
procedure involving use of experimentally determined jack in 
fluences was successfully applied to correction of nozzle profiles 
in a 12-in. supersonic wind tunnel. Considerations of further 


application of this method are discussed 


SYMBOLS 


= longitudinal coordinate along nozzle centerline, meas 

ured upstream from end of nozzle 

coordinate along surface of nozzle, measured up 
stream from end of nozzle 

any finite set of points on the nozzle centerline at 
which an ideal correction is desired 

Mach Number 

flow angle with respect to centerline 

the ideal correction to the centerline airflow which 
would yield perfectly uniform flow 

effect of unit deflection of jack (pair) ion Mach Num- 
ber or angle of the airflow of point x 

symmetric or antisymmetric deflection of jack (pait 
i, positive for deflection of top plate away from axis 
of nozzle 

number of jacks used in the correction 

total pressure of airstream 

stagnation pressure behind a normal shock wave as 


measured by a total head tube 
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INTRODUCTION 


A NOZZLE DESIGNED TO PRODUCE a uniform super 
sonic flow field never produces a flow of perfect 
uniformity because of combined effects of viscosity, 
defects in fabrication, and possible distortions of the 
nozzle caused by air loads. The question of whether 
or not the contour of such a nozzle must be altered to 
improve the flow will depend upon the magnitude of 
the initial nonuniformity, the use to which the nozzle is 
to be put, and the amount of effort required to modify 
the contour. For nozzles used in wind-tunnel testing 
the flow uniformity is usually extremely important, 
particularly if stability and control testing is involved. 

Many wind-tunnel nozzles utilize the solid-block type 
of construction. Modification of the contours of such 
nozzles is accomplished by the addition and removal 
of nozzle material. Various methods of calculating the 
exact distribution of the required corrections along the 
nozzle surface have been given.! In each the same 
basic procedure is followed—1.e., determining the re 
quired flow corrections by measurement of the flow 
in the uncorrected nozzle and relating these to the 
necessary changes in the wall shape by assuming dis 
turbances to be propagated along characteristic lines 
of the flow. The ideal correction, in general arbitrary, 
which would produce perfectly uniform flow, could theo 
retically be made by a single application of this process, 
since the addition and removal of nozzle material could 
also be made arbitrary. It has been found, however, 
that corrections so determined are often less than ideal 
in practice and that the process must be applied re 
peatedly to obtain a particular degree of flow uni 
formity.° 

A considerable number of supersonic wind tunnels 
have been constructed in which the nozzle contours 
are formed by bending thin flexible plates between 
parallel sidewalls by means of jack screws hinged to the 
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A procedure for determining the best possible corre flexil 
tion of a flexible plate nozzle based on the influences , tunm 
individual jacks on the airflow has been given } nti 


Puckett.* Because of certain simplifying assumptions 







































































made in the development of this method, it has no; \NAI 
met with complete success, and most flexible plate no; 
zles have been corrected by various methods of sy 
——— cessive approximation which are partly empirical an Meth 
Via y : partly trial and error. Although satisfactory profiles Th 
a Vel| | , have been obtained by such methods, the procedury estab 
is lengthy and time-consuming, especially when applied irflo 
to nozzles utilizing a relatively large number of jacks fi cer 
An investigation of the method of reference 4 has ind; which 
cated the major shortcoming to be the neglect of th that 
effect of the boundary layer on the process of trans the fl 
a : . mission of disturbances from the wall into the main flow ac 
Fic. | rhe flexible plate nozzle of tunnel E-1 $88 ; : i : ick 
A modification of the method of reference 4, based o simph 
the use of experimentally determined jack influences chang 
plates at intervals. An example of such construction has been found capable of producing optimum corre tial ct 
is shown in Fig. 1, which is a view of the nozzle of tions of a given uncorrected contour with a minimum on wa 
Tunnel E-! of the Gas Dynamics Facility at the Arnold expenditure of time and effort. ire sh 
Engineering Development Center. Corrections to the Tunnel E-1, on which the experimental work con turba 
contours of such a flexible wall nozzle are easily effected nected with this investigation was carried out, is a 12-in wall s] 
by small movements of the various jack screws. The intermittent blowdown wind tunnel with a Mach sectiol 
plate behavior is governed by various laws of bending Number range of 1.3 to 5.0. It incorporates a flexible a 
and arbitrary corrections to the nozzle profiles cannot plate (0.20 in. thick) nozzle that is shaped by a series fluenc 
ye made without using an infinite number of jacks. of twelve jacks spaced at 6-in. intervals. The tunne be det 
I le without using finit ber of jack f twelve jacks spaced at | tervals. The tunnel | bed 
The usual procedure is to correct flexible plate nozzles is supplied with air from a 5,200-cu.ft. high-pressure b 
until a point of diminishing return is reached. This tank and exhausts into a 200,000-cu.ft. vacuum sphere same ¢ 
point is considered to be determined by limitations on Stagnation pressures of 1 to 4 atm. are available, but 
the contours imposed by the laws of bending and by de- operation so far has been limited to atmospheric inlet 
flection of the plate between jacks caused by air loads. pressure. The nozzle itself is a modified version of the 
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FLEXIBLE PLATE SUPERSONIC 


fexible plate nozzle used in the 12-in. supersonic wind 
tunnel at the Jet Propulsion Laboratory of the Cali 
ornia Institute of Technology. 

FLEXIBLE 


INFLUENCE METHODS OF 


PLATE NOZZLE CORRECTION 


\NALYSIS OF 


Vethod Proposed by Puckett 


The method proposed by Puckett? is based upon the 
establishment of curves of influence of each jack on the 
irflow in the test region of a nozzle and the performance 
of certain mathematical operations on these curves 
which theoretically will yield an optimum correction to 
that nozzle. The change in slope of the surface of 
the flexible plates caused by unit deflections of each 
iack can be calculated by the theory of bending of a 
simply supported continuous beam, assuming that slope 
changes for small deflections are independent of the ini- 
tial curvature of the plates. Curves of jack influence 
on wall slope calculated for the plates of Tunnel E-1 
ie shown in Fig. 2. By assuming propagation of dis- 
turbances along characteristic lines of the flow, these 
wall slope changes can be converted to changes in test 
section airflow : 

a) The correspondence of centerline points x; in- 
fluenced by wall disturbances at points x,,, Fig. 3, can 
be determined from a characteristics network. 

b) A change in slope at a point x,, will produce the 


same change in flow angle at x), 
AO(x1) = AP(Xy,) (1) 


c) Mach Number changes at x will be produced by 


the slope change at x,, according to 


AM v} 
tL + [(y — 1) /2].12(x) § AT (x1) 


' 


V A/*(x1) — ] 


d) Proper account must also be made for any waves 
reflected at x,,, as well as for the waves that originate 
at those points. 

Since both the bending theory and characteristics 
theory for small changes in the flow predict that the 
influences will be linear and obey the superposition prin- 
ciple, the net effect on centerline airflow of any combina- 
tion of jack deflections is expressed as a linear com- 
bination of all the individual jack influences, 


> Ad A(x) (3) 
i=1 

For simplicity, attention is focused on corrections 
ol the flow on the nozzle centerline, under the assump- 
tion that a betterment of flow there will result in cor- 
responding improvement elsewhere. To simplify the 
problem further, the deflections A; are restricted to sym- 
metric and antisymmetric deflections of each pair of 
Jacks (top and bottom plate) at the same longitudinal 
position on the nozzle. A symmetric deflection of a 
given jack (pair) will introduce identical families of 


characteristics from each plate. At the intersection of 
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AIRFLOW ——> 








Fic. 3 Propagation of wall corrections to centerline of nozzle 


these characteristics of the nozzle centerline Mach 
Number changes exactly twice those of either family 
will be produced, while the airflow deflections will 
exactly cancel; conversely for antisymmetric deflec 
tions. Inthe above summation (3), if A; is a symmetric 
deflection, the influence /;(v) is the effect on Mach 
Number only; if A; is antisymmetric, /;(x) is the effect 
on flow angle. Thus, centerline corrections of Mach 
Number distribution and flow angularity can be made 
independently. 

All the possible corrections that can be imposed on 
the flow in a flexible nozzle by symmetric or antisym- 
metric adjustment of the 1 jacks is represented by the 
series (3). If measurement of the airflow in the un 
corrected nozzle indicates that an arbitrary distribu 
tion of Mach Number or flow angle corrections, C(x), 
would produce a perfectly uniform flow, the best cor 
rection that can be made in practice will be some ap 
proximation to C(x), since series (3) cannot represent 
any arbitrary function. Assuming that this optimum 
correction will be that which minimizes, in the sense 
of least squares, the differences between the actual 
and ideal corrections over some interval in the test 
region of the nozzle, this condition can be represented as 


“T » 


| C(x) — } Al (x) | dx =aminimum (4 
Ja r=1 
Since this minimum must be with respect to any of the 


n possible jack deflections, 


‘ ‘| [C(x) — >» Al (x) |? dx | = 0 
< x Lva i ] 
(k = . (.) 


When the indicated operations are performed, Eqs. 
(5) reduce to the following system of » linear simul 
taneous equations in the jack deflections, A;: 


n 


z. } T(x) I(x) dx A; = 


1 Le 


ah 
} C(x) I(x) dx Gow 1i1,2...,.8) © 

Puckett noted that the integrals forming the coeffi 
cients of the A; would be fairly tedious to obtain and pro 
posed a variation in which not the best least squares 
correction is sought but one that merely matches the 
ideal correction at a finite number of points. By ad 
justment of n jacks, the ideal correction C(x) could be 
imposed on the flow at any m points x;, which condi- 
tion is represented mathematically by another system of 
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fluence. 
simultaneous equations 


n 
>, Ad (x;) = C(x;) ete... @G 
i=1 

Here the coefficients of A; are simply the ordinates of 
the influence curves at the points x;. By selection of 
the x; near peaks of both the influence curves and the 
ideal correction curve, theoretically the maximum cor- 
rections would occur where most needed, and the cor- 
rections at points intermediate to the x; should be 


small. 


Accuracy of the Influence Curves 


Initial attempts to correct the Mach Number dis- 
tribution of the nozzle profiles of Tunnel E-1 by the 
finite point influence method [Eq. (7)] were unsuc- 
cessful, and an extended investigation indicated the 
failure was mainly the result of inaccuracies in the cal- 
culated influence curves. A comparison of these in- 
fluence curves and similar curves obtained by total head 
measurements in the tunnel at a nominal JJ = 3.5 is 
given in Fig. 4. The change in Mach Number is repre- 
sented by the change in the easily measured pressure 
ratio, P)’/Po, which is a function of local Mach Number 
alone. Since the calculated influence curves would 
introduce errors in the nozzle correction at any point 
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equal to the sum of the errors of the individual influeng, 
curves at that point, discrepancies shown in Fig, 4 cop 
stitute ample reason for the initial failure of the metho 
The discrepancies were initially thought to resy) 
from simplifying assumptions made in applying ¢h 
beam theory. However, measurements of the acty,| 
plate bending characteristics demonstrated that th, 
beam theory was quite accurate. This can be see 
in Fig. 5, which presents a comparison of a calculated 
slope change curve and a corresponding curve give 
measured 


by numerical differentiation of plate d 


flection data. On the influence curves for jack 

(Fig. 4), the discrepancies caused by errors in the bean 
theory can be observed. These are small compared 
with the total discrepancies, indicating that the dis 
crepancies are largely due to deficiencies of the aer 
dynamic theory. Such deficiencies obviously result 
from neglect of the boundary layer. 
boundary-layer effect would be a general upstrea: 


shift of the characteristics network due to both the 


The expected 


curvature of the characteristics within the boundar 
layer and the larger Mach angles accompanying th 
decrease in Mach Experimental influence 
curves for jacks 10, 11, and 12 (Fig. 4) clearly exhibit 
Curves for jacks farther upstream ar 


Number. 


such a shift. 
apparently further complicated by disturbance waves 
originating so far upstream that they negotiate th 
boundary layer more than once in their transit to th 
test section. 

A further discrepancy was noted in the reduction 6 
the peaks of the influence curves. The changes in th 
flow produced at the flexible plate were not transmitted 
in full strength to the main flow because of some damp 
ing effect of the boundary layer entirely distinct from 
the shifting effect just discussed. In most practical 
applications such damping is of relatively minor in 
portance and, consequently, it has been given Itt] 
consideration in currently available treatments of th 
interaction of boundary layers and disturbance waves 
in supersonic flow. Various investigators have demoi 
strated theoretically that a weak wave incident on 
supersonic shear layer does not reflect as a single wa\ 
but as a series of waves constituting a diffusion of the 
wave in the downstream direction. These reflected 
waves, nevertheless, have a net strength equal to th 
strength of the single wave that would be reflected 
in the absence of the shear layer—i.e., the inviscid re 
flection.»-7 Similar results have been found for theo 
retical models of a subsonic layer adjacent to a super 
sonic layer except that the diffusion of the disturbance 
wave extends upstream as well as downstream.” 

Experimental studies of boundary-layer shock-wa\' 
interactions have indicated that at points removel 
from the immediate vicinity of the interaction th 
flow is similar to that which would be expected in the 
absence of the boundary layer. The extent of the 
interaction region, in which severe nonlinearities such 





as separation and reattachment occur, decreases as Ue 
wave strength decreases. It has been observed, 10 | 
example, that the nonlinear region is completely ab 
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FLEXIBLE PLATE 
sent, and the flow is close to the theoretical models for 

wave incident on a turbulent boundary layer at 
f = 3.0. 
cations are that extremely weak disturbance waves 


17 
if Thus, theoretical and experimental indi- 
ze merely diffused and not altered in overall strength 
yy a supersonic shear layer. Reflections that occur 
entirely within the supersonic portion of the boundary 
laver will be free of damping, but waves resulting from 
disturbances at the physical flow boundary will defi- 
nitely be subject to damping between the wall and 
the sonic streamline because of the basic nature of an 
dliptic (subsonic) flow field. This damping would 
be expected to be significant only when the wall de- 
flections are so small as to be of the same order of mag- 
nitude as the subsonic portion of the boundary layer, 
and this provides a possible explanation of the damp- 
ing of the peaks of the experimental influence curves. 
Boundary-layer measurements have not yet been made 
in Tunnel E-1, but at a nominal J/J = 3.50 and atmos- 
pheric supply pressure, the boundary layer is known 
to be turbulent and about 1 in. thick at the end of the 
nozzle. The subsonic thickness of a turbulent bound- 
ary layer at J = 3.50 is of the order of 1 per cent of 
the total thickness, and this is of the same order of 
magnitude as the deflections of the flexible plate used 
in establishing the influence curves of Fig. 4. The 
damping effect at .1J = 2.50 was observed to be less 
than at 3.50. 

The inaccuracies of the calculated influence curves 
could be reduced somewhat by modifying the calcula- 
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puted displacement thicknesses would enable an esti- 
mate of the Mach Number reduction and, in turn, the 
increase in Mach angles. An allowance for curvature 
of the characteristics could be made by use of the ‘‘re- 
flection thickness’? concept suggested by Tucker.® 
The theory of interaction between boundary layers 
and weak waves is not yet sufficiently developed to 
predict accurately the damping effect that was observed 
to occur when the boundary layer is much thicker than 
disturbances causing the waves. In view of the rather 
large amount of work which would be required to pro 
duce only partial reductions of the inaccuracies of the 
calculated influence curves, the readily obtainable ex- 
perimental influence curves are much to be preferred for 
use in making nozzle corrections. 


Superposition and Linearity of Influence Curves 


The series of Eq. (3) can be justified only if the jack 
influences obey the superposition principle and are 
linear functions of the jack deflections. Linearity is 
required from the standpoint of obtaining solutions of 
the correction Eqs. (6) and (7). Simple beam theory 
and characteristics theory for small changes in a non 
viscous supersonic flow both predict that linearity and 
superposition will prevail, but since the boundary 
layer and any interaction with disturbance waves may 
behave in nonlinear fashion for changes of the order of 
magnitude for which the main flow will be linear, the 
extent to which these two requirements are met in any 
























































tions to allow for the boundary-layer effects. Com- particular case must be determined by experiment. 
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P) measurements for a given 
Py values pre- 


A large number of P,’ 
nozzle profile were compared with P,’ 
dicted from some other profile by various combinations 
of the experimental influence curves. Since the dis- 
crepancies between these two sets of data were found 
to represent solely the statistical scatter of the smal] 
experimental errors of the individual influence curves, 
it was concluded that superposition was obtained. 

By direct measurement, linearity of influences was 
found to occur in Tunnel E-1 over an average range of 
+().035 in. jack deflection at a nominal Mach Number 
of 3.5. Departure from linearity was gradual, so that 
for some distance beyond this range the correction 
equations could be solved as though linear and the re- 
sulting jack corrections modified by small increments 
to allow for the The simultaneous 
existence of both linearity and damping of the influences 
is a confusing phenomenon that could be clarified only 


nonlinearity. 


by a more detailed investigation. 


Considerations of Form of Influence Method Most Suitable 

for Tunnel E-1 Nozzle Corrections 

Since it was suspected that the appropriate form of 
influence method for any nozzle correction would de- 
pend upon factors peculiar to a particular installation, 
such as boundary layer and relative jack spacing, a 
number of hypothetical corrections were made to as- 
sumed uncorrected pressure distributions utilizing the 
experimental influence curves of Fig. 4. In Fig. 6 are 
shown several of the indicated corrected pressure dis- 
tributions, differing from those that would be produced 
in the tunnel by only the accumulated random uncer- 
tainties of the individual influence curves. Curve (A) 
is the pressure distribution of an assumed uncorrected 


profile. Curve (B) demonstrates that the finite point 


method, Eq. (7), will produce desired corrections at 
selected points, v;, but overcorrections at intermediat 
points may give a poor overall correction. Curves 
(C) and (D) show the integrated method is a decided 
improvement which increases with the number of jacks 
used. 

The solution for curve (D), however, would have no 
practical value because the correction to the last jack 
was so large that a decided bump would have been pro- 
duced at the end of the nozzle and at the same time 
excessive curvature would have been imposed on the 
plate. Jack that physically _ rea- 
sonable could be obtained by correcting to some pres- 
sure level other than the mean of the uncorrected dis- 
tribution, but such corrections invariably produced 


corrections were 


pressure distributions that were poor, curve (E). It 
was thought that the physical restriction of small cor- 
rections to the last jack before the end of the plate could 
not be recognized by the correction equations, and 
hence it was arbitrarily decided to omit the last jack 
from the corrections. It was later found in practice 
that this omission was unnecessary at low Mach Num 
bers. 

It was also determined that the overcorrections 0! 
the finite point method depend upon the location of the 
x;, being largely eliminated if the «; are chosen to com 
cide with points at which an integrated solution (for 
the same number of jacks) produces the ideal correc 
tion, C(x), exactly. Since these latter points were not 
observed to be correlated in any particular fashion, 
successful x; can be found by trial-and-error process 
only. Thus the integrated method is preferable for cor 
rections involving a small number of jacks (3 to 6 
When a large number of jacks (10 to 20) is used this 


preference may be reversed for two reasons: Furst, 
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becomes inconveniently large (m* + 3n)/2. Second, 


ie number of integrations for the coefficients of (6) 


sensitivity of the overcorrections to the location of the 
will probably decrease because of the decreased rela- 
tive spacing of the points x; at which desired correc- 


tions can be imposed. 


EXPERIMENTAL VERIFICATION OF NOZZLE CORRECTIONS 
UsInG INFLUENCE METHOD 


Procedure 


Attempts were made to correct the Mach Number 
distributions of the nozzle contours of Tunnel E-1 
using the influence method of reference 4 modified in 
accordance with the considerations discussed in the 
preceding paragraphs. Experimental influence curves 
were formed by obtaining the influences for the largest 
jack deflections possible within the linearity limits. 
The percentage error of the influence curves can be 
reduced by this procedure since the probable error in 
total head measurements is the same regardless of the 
size of the jack deflection. For reasons concerned 
with the ease of solution of the correction equations, 
the number of jacks used was restricted to four, and this 
restriction required that the integrated method of 
Eq. (6) be used. Jacks 8, 9, 10, and 11 were those 
used, because jacks upstream of a nozzle inflection 
point were found to have less influence on the detailed 
flow distribution in the test rhombus than jacks down- 
stream and because of the omission of jack 12 for the 


reason discussed previously. 


Po 
Pe 








Example of the Formation and Solution of the Correction 
Equations 


In accordance with Eq. (6), the elements of the fol 
lowing matrix equation were formed by numerical inte 


gration of all the possible products of two experimental 


influence curves for JJ = 3.5, Fig. 4, and C(x) was 
taken to represent correction to Py’ Py = 0.2380: 
34.50 9.¢2 —22.40 —16.07 7S.S7 

2.¢2 24.52 —10.47 —14.04 64.08 to 
—22.40 —10.47 33.82 —6.66] = 16.70 
—16.07 —14.04 -—6.66 37.68 82.20 


All the elements have been multiplied by an appropriate 
constant to give numbers of convenient magnitude 
The solving of Eq. (8S) by the method of minors is 
simplified by the fact that the left-hand matrix is 
always symmetric. The solution yielded the following 
symmetric jack corrections for the top and bottom 
plates, rounded off to the nearest 0.001 in. correspond 


ing to the least count of the jack indicators: 


As = +0.047 in. 
Ay = +0.041 in. 
Aw = +0.034 in. 


An = +0 020 in. 
Pressure and Mach Number Distribution of the Corrected 
Contours 
The above jack corrections were applied to the un 
corrected contour for JJ = 3.5, and the total head sur 
vey then taken in the tunnel exhibited marked improve 
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ment over that of the uncorrected profile, Fig. 7. The 
overall degree of correction was good, but small dis- 
crepancies remained between measured and predicted 
values of the corrected distribution, indicating that 
small errors in the experimental influence curves were 
not completely eliminated. <A similar correction pro- 
cedure was carried out at a nominal JJ = 2.5, and the 
resulting pressure distribution curves are given in Fig. 
S. The influence curves at this Mach Number were 
not obtained as precisely as at 1/7 = 3.5, and this is re- 
flected in a slightly lesser degree of correction of the 
Mach Number distribution. Each of these corrections 
required the taking of only six total head surveys 
i.e., one for the uncorrected contour, four for the in- 
fluence curves, and one for the corrected contour. 

An evaluation of the merit of the influence method 
should be based on the following considerations: 

(1) At JJ = 3.5, the influence method corrected the 
Mach Number distribution from an original variation 
of +0.082 in J/ to a variation of +0.013 in J7 (+0.4 
per cent) utilizing only six total head surveys. A pre- 
vious correction of this contour by a trial-and-error 
process required over 60 surveys to produce a variation 
of +0.014 in J. 
average number of surveys required for nozzle correc- 


This was an extreme case, but the 


tion by the trial-and-error process was 20 to 30 for each 
Mach Number. 

(2) At J = 2.5 the influence method corrected the 
Mach Number distribution from a variation of £0.033 
in ./ to a variation of +0.018 in Af (+0.7 per cent) 
using six surveys. No direct comparison with the 
trial-and-error method is possible because no attempt 
had been made to correct the uncorrected contour for 
this Mach Number by trial-and-error means. 


CONCLUSIONS 


(1) An influence method has been verified in which 
empirical influence curves can be used to determine 
the corrections to the calculated contours of a flexible 
wall nozzle which will produce a more uniform flow. 
Calculated influence curves based on beam _ theory 
and simple characteristics theory will not be sutticiently 
accurate if a relatively thick boundary layer is present. 
Experimental influence curves are so easily obtained 
that they are to be preferred even in cases where the 
boundary layer is thin enough to allow an approximate 
correction for boundary-layer effects. 

(2) Since the superposition and linearity characteris- 
tics will vary from nozzle to nozzle and with Reynolds 
Number and Mach Number for a given nozzle, they 
should be determined experimentally in any given 
case. Jack corrections given by solutions of the cor- 
rection equations should fall within the limits of super- 
position and linearity. To achieve physically reason- 
able jack corrections it was found necessary in some 
instances to omit the last downstream jack from the 
correction. 

(3) Two methods of correcting flexible nozzles by 
influence methods are available: an integral method 
that yields an optimum correction, and a finite-point 
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method that yields exact corrections at a finite numbe, 
of points. The integral method is generally preferah} 
It is, in fact, required for corrections involving a smal} 
number of jacks if reversion to a lengthy trial-and-erro; 
process of selecting the locations of the finite points js 
to be avoided, because the exact location of these points 
is a critical factor in determining the size of overcor 
rections at intermediate points. 

(4) The finite-point method will give an approx; 
mation to the optimum correction of the integral 
method. 
with the number of jacks involved, since the sensitivity 


The degree of approximation will increas 


of the overcorrections to the location of the chose 
correction points will decrease as the relative jac} 
spacing decreases. For nozzle corrections involving q 
large number of jacks, the finite method then may be 
come more practical, since the number of integrations 
required for the integral method increases rapidly with 
the number of jacks. 

(5) The above conclusions have been confined to 
corrections of centerline Mach Number distribution 
by symmetric deflections of the jacks. It is probabk 
that the same conclusions will hold for independent 
corrections of flow angularity by antisymmetric jack 


deflections. 
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Unsteady Aerodynamic Reactions on 
Airtoils in Cascade’ 


F. SISTO# 


Wright Aeronautical Division, Curtiss-Wright Corporation 


SUMMARY 


For one blade in a cascade of fluttering blades, under condi 
tions of small incidence and small camber, the aerodynamic 
e calculated and are found to depend strongly on the 
As compared 


reactions al 
hase between the motions of adjacent blades 
to one isolated blade the lift due to bending and the moment 
jue to pitching are both changed mainly in the out-of-phase 
omponent, whereas the lift due to pitching and the moment 
jue to bending are changed predominantly in the in-phase com 
ponent. These effects increase with decreased blade spacing 
It is concluded that only for exceedingly small gap/chord ratios 
in the effect of cascading account for any large decrease in criti 


| flutter speeds 


NOTATION 


fraction of semichord from midchord to rotation point 
(elastic axis), positive aft 
= blade semichord 
( lift coefficient, C, = L/4q6 


moment coefficient, Cy, = 1/7 /4qb? 


Cos a 
= lift, positive down 
lf = aerodynamic moment, positive stalling, taken about the 
rotation point 
p = pressure 
Ap = pressure difference 
= dynamic pressure, g = pl?/2 
= time 
|" = stream velocity 
= density 
o = velocity potential 
Subscripts 
= airfoil (blade) 
h = due to pure bending 
} = nth blade from central (or reference ) blade 
= wake 
a = due to pure pitch 
Superscripts 


, bar = complex amplitude of a harmonically varying 


quantity 


INTRODUCTION 


i beg POSSIBILITY OF serious aerodynamic excitation 

in the blading of axial-flow compressors first re- 
ceived consideration around 1945. At that time it was 
noted! that certain critical speeds (stresswise) of com- 
pressor test rigs could not be correlated satisfactorily on 
a conventional Campbell diagram: the 
occurring away from the integral engine-order lines. 


resonances 


Received July 9, 1954 
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Inquiries into a possible mechanism have yielded the 
following possibilities : 

(1) Stall flutter, or single-degree-of-freedom excita 
tion due to periodic separation and reattachment of 
the flow about the blades.*~* 

(2) Propagating stall, or the periodic entrance of the 
blades into stalled regions that rotate at a fraction of 
compressor speed around the blade annulus.* 

(3) Cascade flutter due to alteration of aerodynamic 
reactions by interference of neighboring blades.* "|! 

1) Shock flutter, akin to compressibility flutter of 
aircraft components, but complicated by multiplane 


effects and by involved flow passage geometry. ' 
(5) Buffeting between adjacent blade 

to fixed obstructions in the flow passage.‘ 
(6) Mixed effects. 


Calculations of classical flutter for airfoils with com 


rows, or due 


13 


pressor blade properties invariably yielded critical air 
speeds several times larger’ than the critical speed 
observed experimentally. Consequently greater cre 
dence was given to items (1) and (2). However, the 
possibility of strong cascading effects that suitably 
modify the classical flutter expressions so as to account, 
in part at least, for the lowered critical air speeds was 
not overlooked. Mendelson* admitted this possibility, 
and an approximate correction was developed by Belle 
Re 


and 


not and d’Epinay'! for steady cascading effects. 
ports due to Reissner,” Timman,' Gordon," 
Lilley’? may be interpreted as giving the modification of 
classical theory which is dependent on time-varying 
interference effects for the special case of adjacent 
blades oscillating in antiphase. (Actually this par 
ticular phasing corresponds to the wind tunnel wall 
correction problem.) In the present investigation the 
scope is limited to a generalization of these cascade in 


terference effects. 
CASCADING EFFECTS 


With the introduction of cascaded airfoils the flutter 
problem assumes a new complexity characterized by an 
additional length ratio. The dimensionless ratio that 


fixes the cascades geometry is taken to be 


where b} is blade semichord, s is blade spacing, and @ is 
cascade stagger angle (see Fig. | Actually, when 8 
is different from zero, P is complex and contains two new 


length ratios. 
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In addition to this formal complication, certain other 
difficulties of a practical nature appear. In order to 
obtain simple analytical expressions for the lift and 
moment on a blade in a fluttering cascade, it is con- 
venient to assume some regularity or ordered nature 
in the phenomenon. Consequently, the assumption 
is made that the blades are oscillating in the same mode, 
at the same frequency, and with some constant phase, 
mr, between the motions of adjacent blades. 

Restricting the development to a cascade of thin, 
slightly cambered blades and assuming small dis 
placements, small stream turning and plane wakes, the 
integral equation for the amplitude of the vorticity 
distribution, 7,, may be formulated along lines laid 
down by Schwarz,'* Séhngen,”” 7! and Reissner. °* °° 
The result (assuming harmonic motion of circular 


frequency w and a semichord of unity) is 


(jx} S 1), i,(x) = 
| 
~“— Va(é) A[m, P(x — &)] dé — 
oT « 
s fi : a 
rae Yu(é)A[m, P(x — &)] dé (1) 
where the kernel, A, is derived in Appendix (A) and 
reduces, for 8 = O, to 
m| S14), ee oe 
» 


cosh [(1 — |m!)P(x — &)] esech [P(x — €)] (2) 


The vertical velocity distribution at the central blade, 

, is assumed to be given. 

For infinite blade spacing, P > O, Eq. (1) reduces to 
the well-known equation for a single airfoil in which 
K = (x — §)-". 

With the auxiliary condition 


VYw(x) = thle 


where 


*] 
r = | Va(x) dx (4) 
Per 
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k wh / | 
the integral Eq. (1) may be solved for y, given a 
m, P, and @,. Having obtained 7, it 1s then possib) 


to obtain the distribution of pressure difference fro; 
the generalized Bernoulli equation for unsteady floy 


pb p = (Og Ot) + IV(0d Ox 


which, under the previous assumption, becomes 
Aba/pV = Fax) + 1k Valé) dé { 
l 


The oscillatory lift and moment may then be obtaine 
by suitable integrations along the chord. 
APPROXIMATE SOLUTION 


In the present calculations, since A is real and od 
(8 = Q), the kernel is approximated by an expression 0 
the form 
K[m, P(x — &)] = [1/(~ — &)] + 

A(z — & + Bix — t i 


where the fitting constants A and B are determined in 











REDUCED FREQUENCY 





Fic.2. Lift coefficient due to pure bending, Cra, for zero stags 
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Moment coefficient due to pure pitch, CMa, for zero 
stagger, unit gap/chord, and a = —0.35 


Fic. 5. 
) 


the least square sense on the interval (—2P, 2P) for 
each combination of P and m. This approximation is 
With this sim- 
plification the solution proceeds somewhat along the 


integrated only over finite intervals. 


lines of reference 15. Use is made of the Schwarz 


4 


inversion formulas,'* and the lift and moment are ob- 
tained, finally, in the form of three simultaneous equa- 
tions in three unknowns. These equations may then 
(The third unknown is 
The 
development is lengthy'’ and is omitted here. For 


be solved by Cramer’s rule. 
the second moment of the pressure difference.) 


reference purposes the results may be displayed in the 
following contracted form: 


L = rpbV?A- (hha + hoa& — hoh/d) (S) 


VW = wpb? V2A (ne + joe — joh b) — abL (9) 


where L, 7, a, and h are the complex amplitudes of 
lift, moment, pitching displacement, and bending dis- 
placement, respectively. A, 4, Ao, ji, and je are com- 
plicated algebraic functions of A, B, k, and C,. The 
“cascade Theodorsen function,” C,, is given by the ex- 
pression 


Ce. = (¥1 + fids)/(V1 + fda + i(Yo + fiJo)] (10) 


where the Y’s and J’s are Bessel functions with argu- 
ment k. A listing of the functions A, f,, 4, /e, 7;, and 
j2 appears in reference 18. 
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This 


asymptotic behavior as blade spacing increases (P + 


latter group of functions shows the proper 


Q). Furthermore, a comparison with an exact soly 


tion (for m = 1, P = 1.82, reference 17) shows th 
accuracy of the approximation, Eq. (7), to be fairh 


good. 


DISCUSSION OF RESULTS 


Expressions (8S) and (9) are rather cumbersome whe: 
expanded. Consequently, the discussion will be cor 
ducted for some numerical examples in which 
— 0.39. 

The first numerical calculations are for a cascade oj 
zero stagger and pitch/chord ratio of unity. The vari 
ous aerodynamic coefficients are displayed graphically 
in Figs. 2 through 5 as functions of phase mz and r 
duced frequency, k. In each case the curve for ar 
isolated blade, if it were drawn in, would fall approxi 
mately midway between the curves for m = 0 and m = 
(with some minor exceptions). 

It is interesting to note that in no case is there a 
radically different type of behavior due to cascading 
The lift due to bending and the moment due to pitching 
(C,, and Cy,,,) are both changed mainly in the out-of- 
phase component, whereas the lift due to pitching and 
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Fic. 6. Lift coefficient due to pure bending, CLa, for zero st4e 
ger, several phase angles, and several gap/chord ratios. 
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the moment due to bending (C;,, and Cyy,,) are changed 
predominantly in the in-phase component. 

The effect of solidity, s/2b, on the oscillatory aero 
dyvnamie reactions is illustrated by one sample plot, 
Fig. 6, showing lift due to bending. Only for the ex- 
eedingly small gap/chord ratios do the aerodynamic 
reactions seem to be changed sufficiently to account for 
mv several-fold changes in critical flutter speed. This 
statement is confirmed, essentially, by the experimental 
work of Lilley.!” 

It seems safe to conclude, therefore, that the radically 
lowered critical speeds actually encountered in axial- 
fow compressors are not caused by interference effects 


in classical flutter aerodynamics. 


APPENDIX (A) 


In this appendix an expression is developed for the 
velocities induced by a cascade of fluttering blades 
under conditions of two-dimensional potential flow. 

Consider a row of equidistant point vortices, of spac- 
ing s, and of strength 7,,(, ¢) d&, disposed along a line 
inclined at an angle 8 to the y-axis as in Fig. (7). The 
velocities induced by this array at any point z may be 
expressed, according to the theorem of Biot and Savart, 


in complex notation, as 


: , 1 dé vné, 1) 
dw(s) = du —idz : 
23 z-— s, 


where 
So = § + inse””, = @, £1, +2,... 
ind zis the complex variable, 2 = x + ty. All summa- 


tions in this appendix are on n, and the range of this 
index is from — to +o, unless otherwise noted. 


If interest is confined to velocities at the x-axis, there 


° 38 + + 
te” dé yn(é, t) 
dW\xX) = 
2s mre'’(x — ¢)/s —t 
“as TE AX s . int 


In order to synthesize a cascade of blades it becomes 
necessary to distribute vorticity along lines representing 
the chords and wakes of the airfoils. The expression 
lor the total induced velocity at any point on the x-axis 


4 2. yn(é, t) dé 
iy = 
27 1 P(x — t) — mr 


s and the semichord 6 has been taken 


becomes 

w(x) =u — 
where P = rhe’ 
equal to unity. 

Assuming now that the (distributed) vorticity is an 
harmonic function of time, say 

ynlé&, t) = yo(&) cos (wt + nm), (jm! = 1) 

the following summation appears under the integral 
sign 


c=) 


+ . ( 
dX (cos (wt + nmm)/|P(x — &) — inn); 
rhis summation may be expanded, in several steps, to 
vield 


ACTIONS 


( 
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f l 
o = COs al \p — ie py cos numa X 
P(x —é 
1 {i a 
; : Pe sin wl - r B sin mn x 
P(x — &) — inn inn 4 { 
| > sin m1) 
P(x — &) ine inn | inn 4 


where >>’ indicates that m = 0 is to be excluded in the 
summation, and in the first summation the terins for 
n and —n are always to be taken in pairs. By an ap 
plication of the Mittag-Leffler theorem,** it may be 


shown that 


= l l 
> COS NM 1 = = = 
P(x — &) — inn int 
m ' 
cosh —m) P(x é x 
m 
h [P . 
cscn [F(x — &) 
be Pix — t 
= l l 
S°’ sin nmr : : + - = 
P(x — t) — intr int _} 
—_ m . ; 
1 sinh —m) P(x — £)| X 
m 
: _f[\im 
esch [P(x — &)] — 1 m 
m 
and consequently 
m , : 
ao = cos wt cosh —m) P(x —&)| X 
m 
esch [P(x — &)] —7 sin wt X 
m : | ; : 
sinh —m) P(x — &)}esch [P(x — &~ 
| \ m J 


since it can be verified that 


>’ sin mmx/ine = —i[( mm m| 














Fic.7 Nomenclature for Appendix (A 
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The final expression for the vertical component of the 
induced velocity may be written 
yo(E) X 
[R.P.(¢’) cos wt + 1.P.(0”) sin wt] dé 


where 


m 
o’ = P cosh ( —m]) P(x — &) | esch [P(x — &)] 
m 


m 
—m) P(x — &) 
m 


For P real and with the replacement of yo cos wt by 


o” = P'sinh esch [P(x — &)] 


je’, this reduces to Eqs. (1) and (2). 
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SUMMARY 


The theory of the steady flow of a viscous compressible fluid 


i flat pl it¢ 
8 js extended by a more complete discussion of the flow in 


it high Mach Number due to Lees and Prob- 


inviscid laver between the shock wave and the boundary layer. 


shown that similar solutions exist in this laver, analogously 


those found by Li and Nagamatsu?® in the boundary laver, and 


t the two may be joined to give, allowing one minor assump 


ta full account of the flow. It is shown that the boundary- 


r equations may be reduced to those for an incompressible 


uid and that the von Karman-Pohlhausen method describes 


flow in it with good accuracy The tangent wedge approxi- 


tion for the pressure on the plate, used by Lees and his col! 


orators, is found to be in deficit by 10 per cent for air. Finally, 


s shown that the theory for weak interaction’ cannot be ex- 


led further without a complete knowledge of the flow. 


NOTATION 


distances along and normal to the plate 
= velocity components in the (x, y) directions 


pb, p, 1 = pressure, density, and absolute temperature at 
YY 
I itio of specific heats 
Prandt] Number taken = 1 
= kinematic viscosit\ 
velocity of undisturbed fluid 


= velocity of sound in undisturbed fluid 


= | 
y=} equation of shock wave 
y stream function 
S(y entropy function 


Subscripts 


= conditions in undisturbed fluid 


= conditions just behind the shock wave 


= conditions at the outer edge of the boundary laver 





(1) INTRODUCTION 


, A PREVIOUS PAPER,!? which we shall refer to as I 
throughout, the flow induced by the impulsive 
motion, in its own plane, of an infinite plate immersed 


Rect iv ed April 23, 
*Rese arch 


England 


1954. 


Fellow, on leave from The University, Bristol, 


‘In a paper to be published in this JouRNAL, Mr. Lester Lees 
onsiders the effect of this assumption and concludes that the 


O[( Aly 


expect from the unsteady problem, and > 0 as I —> 


ire 


“ror incurred is indeed (Ux)] /* Y as we would 


©, keeping 





Ux finite 10 and 


in comparisons with experiments at this Mach 


He also estimates the error when JJ = 
concludes that, 


Number, it cannot be neglected. 


On the Motion of a Flat Plate at High Speed 
ina Viscous Compressible Fluid—II. 
Steady Motion 


K. STEWARTSON* 


Guggenheim Aeronautical Laboratory, California Institute of Technology 


in a viscous compressible fluid was considered, assum 
ing that the characteristic Mach Number .\/ was large. 
It wa 
where vy 
the velocity of sound in the undisturbed fluid, .1/ 
U doy, 


and ¢ is the time measured from the instant at which 


J 
s shown that, provided [yp (ao7t)] log J < 1 


is the coefficient of kinematic viscosity, do is 


where U is the constant velocity cf the plate, 
motion began, it is possible to justify, on the grounds 
of consistency, a division of the flow field into three 
extending from a shock wave 


regions. In region (a 


outwards, the fluid is undisturbed. In region (c) ex 
tending from the plate part of the way to the shock 
wave, viscous effects are important, but the pressure 
changes across it are small so that boundary-layer 
theory may be used to describe the flow within it. In 
this region the gas is extremely hot in comparison with 
the fluid in (a) but has negligible density. The outer 
bound of (c) is clearly defined, and between (a) and 
(c) there is a third region (b) in which the effect of vis 
In large 


entropy gradients occur near the junction with (c 


cosity is negligible. this region, however, 

A similar situation has been shown to hold for the 
steady motion at high Mach Number past a semi 
infinite flatplate.. On atypical streamline coming from 
infinity, the flow field is undisturbed until it meets a 
shock wave emanating from the leading edge of the 
plate. Thereafter the streamline is distorted under 
the laws of inviscid flow, but at no stage do local vis 
cous affect the fluid Between nearly 


all the streamlines and the plate is a boundary layer 


forces on it. 
containing an infinitesimal amount of fluid, but, be- 
cause of the enormous temperature in it, it is of con 
siderable size and of importance in determining the in 
viscid flow. 

Two parts of the flow are amenable to treatment, and 
both have been considered by Lees and Probstein.® ° 
First, near the leading edge, but sufficiently far from it 
for the boundary-layer approximations in the viscous 
layer to be valid, there is a strong interaction field in 
which the shock wave exerts a large effect on the flow 
through it. Second, a long way downstream there is a 
weak interaction field in which the shock wave exerts 
only a small effect on the flow through it, so that be- 
hind it there is only a small perturbation of the uniform 
conditions upstream of the shock. The region in which 
the shock wave is strong was first considered by Li 
and Nagamatsu,*® assuming that the shock wave coin- 
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cided with the outer edge of the boundary layer. 
Their work is of importance because they demon- 
strated the existence of “‘similar’’ solutions in this re- 
gion, although, as Lees® has pointed out, the basic as- 
sumption leads to a contradiction, because the edge of 
the boundary layer must be a streamline; therefore, 
from the Rankine-Hugoniot relations it cannot also 
be a shock wave. Lees®** improved their theory by 
showing that there must be an inviscid layer between 
the shock wave and the boundary layer, and, giving a 
partial justification of the matching procedure, he 
adopted to join the two layers. He did not discuss 
the inviscid layer in detail beyond pointing out that 
conditions at the edge of the boundary layer are deter- 
mined with good accuracy by means of the tangent 
wedge approximation; nor did he show that the two 
layers were contiguous. 

The main purpose of this paper is to show that there 
is a ‘‘similar’’ solution in the inviscid layer which may 
be joined up to the “‘similar’’ solution in the boundary 
layer without inconsistency so that there is no inter- 
mediate layer. In the course of the argument a slight 
change in the boundary condition at the edge of the 
boundary layer has to be made. Although it has not 
been possible to show that such a change has no appre- 
ciable effect on the flow field in the steady problem,* 
this assumption is made more plausible by the fact that 
in the corresponding unsteady problem (1), a similar 
change can be completely justified. Using our solution 
valid from the stock wave to the plate, it may be de- 
duced that when the shock is strong the pressure on the 
plate is in deficit by 10 per cent if the fluid is air. In 
addition, we show that the von Karman-Pohlhausen 
method may be used to describe the flow throughout 
the entire boundary layer with an error varying from 
—1 per cent to +2!/. per cent. Finally we show that 
the solution obtained by Lees and Probstein’ when the 
shock wave is weak cannot be further improved without 
a knowledge of the entropy variation across the in- 
viscid layer—i.e., without a knowledge of the entire 
inviscid flow. 

There is little doubt that the Navier-Stokes equations 
are inadequate to give a full and correct description of 
the flow field and that, especially near the leading edge, 
molecular considerations will be needed. However, 
it is likely that the assumption that the fluid is a con- 
tinuum will be justified if M[v)/(Ux)]’* <1," where 
u is the velocity of the fluid at infinity and x is distance 
downstream from the leading edge. As we shall see 
below, this is precisely the condition under which the 
assumptions of this paper are justified. Further, 
although the ‘“‘inert’’ degrees of freedom would require 
a longer distance than that to respond to changes in 
conditions, it is probable that in any practical problem 
the temperatures would be sufficiently low to assume 
that either they are not excited at any time, or, if they 
would then have little influence of the flow. Finally, 
it is pointed out that the molecular considerations, as 


* See footnote on page 303. 
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typified by the existence of slip flow, may be inter 
preted as providing some sort of lubrication between th, 
plate and the main stream, reducing the impact of th, 
one on the other, so that the results of the conventiona) 
boundary-layer theory would give upper bounds for th, 


disturbances. 
(2) THE EQUATIONS OF MOTION AND THE Inviscy 
LAYER 


We assume that the plate occupies the positive ha] 
of the x axis and that the fluid is streaming along jt 
The Mach Number 


ado < 1, where dp is the velocity of sound at in 


with uniform velocity u at infinity. 
M=U 
finity and the region in which the flow is disturbe 
is bounded by a shock wave and immediately beloy 
the flow is inviscid. We are not concerned here wit} 
the thickness of the shock wave nor with the nature of 
the flow within it, beyond assuming that it is thi 
and that conditions on either side of it are determined 
by the Rankine-Hugoniot relations. Under the as 
sumptions of a continuum, these properties of the 
shock wave may be justified. We denote by p th 
pressure, p the density, 7 temperature, y distance 
normal to the plate, u, v the components of the fluid 
velocity in the x, y directions, respectively, and y the 
ratio of the specific heats of the fluid. Let the equatior 
of the shock wave be 

Y = V(x) (2.1 
where Y is at present unknown; but we shall assume 
that (dY/dx)*? < 1. If subscripts 0, 1 denote the u- 
disturbed values and the values just behind the shoe 
wave of the variables, we have, following Goldsworthy, 


u— U = UOMdY/dx)? 


o, = [2U/(y + 1)] dY/dx) X 
(1 — [1/M2(dY/dx)2]} 


l 


Pi/ po = [2¥/(yv + 1) Med Y/dx)? — = 
Ky — D/(y+ 


pr/oo = [(y + 1)M2(d¥/dx)*]/[(y — 1)M? X 
(dY/dx)? + 2] 


The flow behind the shock wave is now discussed 
adopting an heuristic argument. We assume that 
the flow is either inviscid or of the boundary-layer typ 
This 
procedure was used in I with success, and a consider- 
able amount of space was used in showing that the pro- 
Here we shall not devote much 


in which the pressure is a function of x only. 


cedure is consistent. 
space to such a justification, merely pointing out the 
For ex 


2.2) are 


circumstance in which the argument is valid. 
ample, the shock conditions given in Eqs. ( 
only valid if dV /dx < 1, which implies that } <.%. 
Consider the following system of equations. 
[pu(Ou/Ox)] + [pv(Ou/Ov)] = —(Op/dx) + 
(0/Oy) [u(du/dy)] (2.3 


pu(dv/Ox) + pv(dv/Ovy) = —Op/Oy (2.3 
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MOTION OF A FLAT 


O Ox) (pu) + (O/Oy) (pv) = O (2.5 

o7 oT Op Op 
1 + MW pt’ —: — 9 = 
“ov |” Oy ox oy 
(o) O (“ ~ OL ) 
Mm + (2.6 
oy Oy \a Oy 
p = RpT ER, 


where c, is the specific heat at constant pressure and 


,is the Prandtl Number. Eq. (2.5) is the equation of 


mtinuity valid everywhere. Eq. (2.5) is from the 
component of the momentum equation valid in the 
undary layer regime and in the inviscid regime. 
Eq. (2.4) is from the y component of the momentum 
uation valid in the inviscid layer. In the boundary 
laver its effect on P must be negligible. We shall show 
this below. Eq. (2.6) is the energy equation valid in 
oth regions, and Eq. (2.7) is the equation of state 
valid everywhere. There are many terms neglected 
in these equations, but it may be shown that in both 
regimes their effect is negligible if ¥ <x. 

Throughout this paper we shall assume that the fluid 


possesses the properties 
og = | 
gp = pelt / To) 


No fluid has these properties, of course; but if we assume 
them, the discussion is considerably simplified, and the 
solution in this idealized case will help us to understand 
the behavior of more realistic fluids. 

The boundary conditions satisfied by the fluid are 
that just behind the shock wave Eq. (2.2) is satisfied 
while on the plate y = 0, x > 0, 


oT /oy = 0 (2.9) 


u=v =), 


taking it to be thermally insulating. The equations 


f motion then have a simple integral—namely, 


— 1)/(2a9?)] (U? — u? — v*) (2.10) 


T/T; = 1+ [(4 


Further, if we define a stream function y by the equa- 





tions 
py(Ow Ov) = pu; pi(OW/Ox) = —pr- (2.11) 
Eqs. (2.4) and (2.5) reduce to 
Ov OX = — (1/ po) (Op/ Oy) —_ 
(O/Ow) (v/u) = (0/0) (po/ pit) _— 


It is noted that y is a continuous function of x and y 
equal to Uy in the undisturbed fluid. Just behind 
| the shock wave, u — U = UO(dY/dx)?, so that u — U 
| may be neglected in comparison with U’ and in Eq. 
*.12) « may be replaced by U’. The relevant equa- 
tions in the inviscid layer behind the shock and also 
the conditions at the shock are similar in the unsteady 
For a 





problem, and the solution is of the same form. 
luller discussion of the relation between the unsteady 
and the steady problems, the reader is referred to a 
Paper by Goldsworthy.* It should be pointed out, 
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however, that there is a misstatement in it, for he says 
that the boundary layer may be neglected in a discus- 
sion of the inviscid field. It is now well known, and 
will be shown here, that such is not the case. The 
similarity between the unsteady and the steady flows 
breaks down either when (# — UL’)/u becomes of order 
one or when the inviscid assumption is no longer 
valid. It may be shown that both of these occur simul- 
taneously. From the energy equation in the inviscid 
regime we find that p, po = S(x) (p/ po)’, where S is a 
function of y only, which may be found from its value 
just behind the shock wave. If we assume that the 
equation of the shock wave is 

Y = Oe (1 + ac +... (2.13 
leaving the justification until the next section, we 


find that 


vote (= (4) 4 (8) 
all S(y +1) \y + 1 Ud? Uv 


and, if we write 


» = ¥/UY 
v = [(2U/(y + 1)] (dY/dx)V(y, x 
p/po = [2yM?/(y + 1)] (dV /dn)*P(q, x 
V = Voy) —*x |; ali(n) + 


16 V | - 
oln) + .., 
gM27g: 7" 


1 
= Py(n = 2 l 


16 
P> } athe 
gieg2 2" | 


the boundary conditions and equations satisfied by the 


aPi(n + 


P’s and I!s are the same as for the same symbols in 
I, Section 3. We 
inviscid solution is valid anywhere behind the shock it 


may conclude therefore that if the 
is valid everywhere in ¥ > 0, because the only singu- 
larities develop there. In particular the proportional 
change in # remains small in y > O, but it has a singu- 
larity at Y = 0. The boundary layer in terms of y can 
therefore only occupy a narrow layer near y = 0. 
Since Op /Oy is bounded as y > 0 + and v increases 
from the plate to the edge of the boundary layer, it 
follows from Eq. (2.11) that the change in ~ across the 
layer is small, and so a conventional boundary-layer 
theory may be used. It is noted that the end of the 
inviscid layer is associated with a rapid increase of 7 
and rapid decrease in uw and p, all of which we must 
have for a boundary-layer theory. Further, since y¥ 
is a measure of the mass of fluid crossing the shock 
wave, almost the whole mass of fluid confined between 
the shock and the plate is in this inviscid layer, and the 


mass in the boundary layer is negligible. This does not 








306 JOURNAL OF THE AERONAUTICAL SCIENCES—MAY, 1955 


mean that the boundary layer is of negligible thickness 
in the physical plane however. In terms of y, the 
boundary layer is nearly as thick as the inviscid layer 
when x is small. This is shown below. In the theory 
of Li and Nagamatsu,’ it was assumed that the shock 
wave coincided with the outer edge of the boundary 
layer, which is equivalent to assuming that the whole 
of the mass behind the shock wave is confined inside 


the boundary layer. 


(3) THe BOUNDARY LAYER 


It was shown in the previous section that once the 
flow just behind the shock wave was admitted to be 
inviscid it followed that the viscous layer was confined 
Further, 


to the immediate neighborhood of y = 0. 
the variation of p across this layer is small, so that in 
it p may be assumed to be a function of v only. The 
relevant equation of motion is then, from Eqs. (2.5) and 


(2-01), 


Ou 1 dp uO [up Ou 
u = — + - (3.1) 
Ox y p dx Po oy Po oy 
with boundary conditions vu = 0 at y = Oandu—> U 


outside the boundary layer. It is noted that this 
boundary condition may be replaced by u —~ U as 
y — ~ because y is still small in comparison with Y 
when wu is not sensibly different from Ul’. Further, 
the boundary condition on u as y > © is really u — 
U~ Uy” *’, from the behavior of 7 in the inviscid layer; 
but it is likely, although we have not proved it, that 
the simpler condition above is justified.* It is cer- 
tainly justified in the unsteady problem to which the 
reader is referred. Using the energy integral in which v 
is neglected in comparison with uw and assuming yp « 7, 
the equation for u reduces to 


Ou 1 po y=]... a» \ dp 
= — l U?— uv’ 
" Ov po Pp ( sl 2ao” ( . ) dx T 


p ra) ( *) 
Vo u u (3.2) 
Po Oy Oy 


We now write u* = u(po/p) , whence 1* satisfies 


(y—1)/24 2 on aeaaes 
( p ) at ou* _ _ (: 42 't) dp 4 
po Ox yp 2 dx 
57 —3)/2) 

sis a) ou* 
Y% ( P| u* (us : ) (3.3) 

Po oy oy 
In this equation we may replace [(y — 1)/2]A/* + 1 


by [(y — 1)/2]12*. Then if we write 


z (By —1)/27 
oS = [ (*) dx 
J0 po (3.4) 


W = U(po/p)~? 
ou* 1W ra) ou* 
u* = = i + wu* oy (u* =| (3.5) 


which is the same as the boundary-layer equation for 


* See footnote on page 303. 


an incompressible fluid. The boundary conditions 


u* = Oat y = Oandu* > Wasyo> 


The normal procedure to determine the z componer 
of the velocity at the edge of the boundary layer js ; 


use the formula 
o (> 


v= — u dy 
. Ox J, . ‘ 


But this method fails here because at the edge of tj 


boundary layer p ~ y 
for all values of x. 
which 
0 (* pody 
Ox JO pu 


v= Uu 


a) im po 
u dy X 
Ox J 0 pu 


Now in the boundary layer | 


uw”) and may therefore be neglected. Hence 


y—1 ra) o ££" sf U i 
y= Meu ; | ( _ .) dy | 
2 OrL p Jo \u l | 


and to find the value of v at the edge of the boundar 


layer we may replace the upper limit by 


y is large on the scale of the boundary layer u — [’; 
At the edge of the boundar 


exponentially small. 
layer therefore 


ae ox 
BCG-a4 
p Jo u Fy v 


since u = U there. 


the plate in the corresponding incompressible boundar 


layer so that 
Oy /oy* = u* 


Eq. (3.8) may be written 


Y- l ca d |(S) 1)/2 
Vz = M?*l 
. 2 dx p " 
we) 7 
J0 WW? a ae 


Further, in the boundary layer 


vy - = | * 3 1)/2 
y= | a 3 ur (2*) x 
J0 pu 2 p 


using the same arguments as before. Hence t 
boundary-layer thickness ly is definite and 
Ve = U(dY «g/dx) oe 


in marked contrast to the incompressible problem whe! 


and is zero on inviscid theor 
Instead we use Eq. (2.12) fro 





©, since wher 
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MOTION OF A FLAT 


the outer edge It is 
of interest to note that according to Eq. (3.10) the 


of the boundary layer is a streamline of the 


to the boundary layer is vague. 


outer edge 
inviscid flow, and this is in agreement with the con- 
clusion of the previous section that the mass of fluid in 
the boundary layer is negligible in comparison with 
that in the inviscid layer. 
If the pressure in the boundary layer is 
b/pb Onx 


P/Po = & 


where Qo is a constant, then s ~ wx and IV ~ 


- so that if we write 


Widf d¢) 


“” = 


, solution of the boundary-layer equations may be 


found in which / is a function of ¢ only. In fact 


les) + [f(d*f/de*)] + 


l(y — 1)/y] [1 — (df/dg)?] = 0 (3.11) 
/ df/dé = Oat¢ = 0 
and 
df/di > 1las¢— 
Eq. (3.11) is one of the family found by Falkner and 


Skan and integrated by Hartree,* from whose paper the 
value of v at the edge of the boundary layer, 


x / Vo 
is a) nus x 


CD -G)]« « 


We note that this value of vp is 


12) 


may be calculated. 
of the right form to effect a match with the inviscid 
solution discussed in the previous section. The sig- 
nificance of the similar solution of the boundary-layer 
equation was first pointed out by Li and Nagamatsu,’ 
and it was obtained in this form by Lees.® The similar 
In general, p may 


of which that 


solution is only valid when p «x 

be represented by a power series in x 
is the leading term. The appropriate solution of the 
boundary-layer equations is then found in the form of 
another series in powers of x *of which the leading term 
is given above. In view of the complicated nature 
of the transformation such a procedure would rapidly 
become unwieldy. Since we would expect the pressure 
on the plate to decrease uniformly as x increases, it 
follows that IV” should increase with z, and therefore 
the von Kérman-Pohlhausen method applied to Eq. 
Following Goldstein? we 


3.0) Should be accurate. 


define 
(3.13) 


vo) (dW /dz) 


A = (8? 


Where 6 is the boundary-layer thickness in the von 


Karm4n-Pohlhausen sense. Then 


PL 
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"(1 — 2) aye - 
Jo We} - 


6(0.4175 — 0.00939A — 0.000110 A? 6bA(A 3.14 
and 
da d ; 
g(A (log 17) + [A*h(A) + A] X 
dz az 
d ( = 5 
log 3.19) 
dz dz 
where g and / are tabulated in reference 2. If now we 
transform back to p/p» and x, we find that 
y-—1 — d 
Ve = AM? (Uv) x 
2 dx 
A“K(A 
3.16 
{(y — 1)/2y] (d/dx) (—p/po)} 
where 
da y¥— 1 37 l 
ee : g(A) + — [A’L(A rab x 
dx { 24 + { 
d 4 d d p y 
log p/po + [A’h(A) + A] — log (3.17) 
dx ax dx Pp 


The initial value of A, Ao is found by setting p_= Qo 


and requiring that A» should be finite. It then follows 
that the right-hand side must vanish whence A 
3.256 if y = 1.4, leading to a value of v at the edge of 
the boundary layer equal to 


0.391 APU [ye (QoU)] * a 3.18) 


If we use the exact boundary-layer solution |Eq. (3.12) ], 


the numerical constant is found to be 0.594 so that there 
is a deficit of about 1.0 percent. The form for p,«x~ “’, 
assumed above is valid when x is small, and, when x is 
large, p is nearly equal to py. In this case the value 
of v at the edge of the boundary layer, calculated by the 
the 
per 


excess ol 


Thus the 


von Karman-Pohlhausen method is in 


correct answer by about 2'/s per cent. 
centage error involved on assuming vg, as given by Eq. 
(3.16) varies from —1.0 per cent to +2'/. per cent, 
and it is plausible to assume that in between the error 
is sufficiently small for us to take Eq. (3.16) as a suth- 
ciently good approximation. 


If the pressure is of the form 


Qo + Q, + O(x + 


p Po = 


we find that when y = 1.4 


ve = M2U[v/(QoU)] x “ [0.391 
0.474 (Q,/Qo)x ° +...) 


(3.19) 


We may now effect a match between the inviscid solu- 
tion at y = 
yy = «. From the previous section and from I, Sec- 
tion 3, we find that if y = 1.4 and y = 0 


0 with the boundary-layer solution at 


(l+x 
[4.748/ (AP li aa 


b/ po = 0.393.M282x 16.984a - 


3.20) 


Ve = 0.444U 3x ‘(J +x” , 3. 19Sa 


(6.218/( M1782) ]; +...) 
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Hence, comparing Eqs. (3.19) and (3.20) we find that 


M4 


3 = 1.188(Al*9/2) 


uo (3.21) 
a= 0.008 | . | 
AL® v9 


whence the pressure on the plate is 


Pw/ po = 0.555 M*[m/(Ux)]” + 0.489 +... (3.22) 
and the equation to the outer edge of the boundary 
layer is 
0.704x [(AL®9)/(Ux)] “4 — 

().460x [(M1%y) /(Ux)]7 "7 +... 


v= Vy 
(3.23) 
It is fairly easy to see that the pressure on the plate is 
a funetion of [(.1/%)/(Ux)]* only. Lees and Prob 
stein’ showed that when x is large the pressure on the 
plate is given by 
1 + 0.349 [(Al®m)/(Ux)] 7? + 
0.052[(.1/%)/(Ux)] +... 


Pu Po = 


(3.24) 


and the two solutions may be joined in an approximate 
manner by making use of the tangent wedge approxi- 
mation in the inviscid region, as suggested by Lees and 
Probstein,* in conjunction with the von Karman-Pohl 
hausen method in the boundary layer. The tangent 
wedge approximation has been used in this way in I 
where it gives only fair agreement with the exact solu- 
tion near x = 0. It is felt accordingly that there 
would be little point in using the method here.  Per- 
haps the only suitable method is the method of char- 
acteristics, already used by Eggers and Syvertson! in 
related problems. 

We can now determine the region of validity of Eqs. 
(2.3) to (2.7) and of the division of the flow field into a 
boundary layer and an inviscid layer. We need not 
discuss it in detail however, because this has already 
been done for the unsteady problem in I. The equa- 
tions of motion are set out by Howarth,°® and we find 
that we are entitled to reduce them to Eqs. (2.5) to 
(2.7) if (Y/x)? < 1—.e., if 
M?(v/Ux) <1 (3.25) 
Further, we may divide the flow field into inviscid and 
viscous fields if, in terms of y, the thickness of the 
viscous layer is much thinner than that of the inviscid 
field. The condition for this is also Eq. (5.25). We 
conclude therefore that the results of this paper are 


valid in the range 
M2(w/U) <x < M*®(n/U) 
and note that, according to Tsien,'! the lower bound 


is also the lower bound for the validity of the assump- 


tion of a continuum. 


(4) THE INFLUENCE OF THE ENTROPY 


In the previous section we obtained the first two 
terms of a series for the pressure p, on the plate in 
. o 1/» . . - 
ascending powers of x “’, and it is clear that any further 
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terms may be obtained if desired. Lees and Prob 
stein’ have bound the first three terms of a. series 
for py , neglecting entropy 
variations in the inviscid layer [Eq. (3.24)]. We ey 
amine here the effect of the entropy variation at larg, 
values of x and show that, notwithstanding the singy 


in descending powers of x 


larity in the entropy at the edge of the viscous layer. 
its effect on Py is O(v 3ut to calculate the fourth 
term of Eq. (3.24), which will therefore depend on th 
entropy, we need to know its distribution in the entir 
inviscid region. It follows that no more terms of th 
asymptotic series for p, can be found until the com. 
plete solution is known. 

To show this we reconsider the equations of motion 


in the inviscid layer—-viz., 


v/ ON —(1/po) (Op/ Oy) 
1.1 
Ov/Oy = —[S v(p/ Po) | (0/Ox) (p Po 
The boundary condition at Y = 0 is that v = v,(x) and 
at the shock wave Eq. (2.2) must be satisfied. Given 


vy, we have three conditions to determine— p, vz, and 


y—and then vy, may be found in terms of p, from the 


boundary layer, giving us an equation for p,. The 
method adopted by Lees and Probstein’ was to neglect 
| — S and use the theory of isentropic flow. They 
found that the entropy at large values of y and x is 
«yp *, because the pressure rise across the shock wave 


<y  *, and concluded that the entropy rise across the 
shock wave would not affect the first three terms of 
Eq. (5.24). There is however the entropy variation 
at small values of y to consider, and we can treat the 
leading terms due to it by replacing the coefficient of 
Op/ox in Eq. (4.1) by —[1/(ypo)] (1 + R), where 
1+R = S', and by neglecting entropy variations 
across the shock wave. Contributions to the fourth 
term of Eq. (3.24) arise (1) from the isentropic solution, 
(2) from the entropy rise across the shock wave, and (8 
from the entropy variation near Y = 0. These three 
contributions are all independent, and the first two 
may be deduced from the asymptotic solution already 
We now calculate the third using the simpli 
The correction that would 


known. 
fied form of Eq. (4.1). 
have to be applied would not affect the fourth term of 
Eq. (3.24). 

We write 


v = ve[x — (¥/ao)] + % 
Pp = pot Ao pode [Xx — (ap ao) | + pi 


where we assume and justify a posteriori that 1 < % 
and py < dyppty. The boundary conditions are that 
%, = Oat y = Oand p,>O0as ym o~-. The leading 
terms of Eq. (4.2) are simply the leading terms of their 
asymptotic expansions. Entropy variations near ¥ = 
0 only affect p and v via the interaction between S and 
these terms. 

Substituting for p and v in the modified form of Eq 
(4.1) and neglecting ),R in comparison with a poe, 


we obtain 
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2p 1 O°py pr ” / 
 .. P = Rip)ve" [x — (p/ao) ] (4.3) 
oy? ay” Ox* ao 
with boundary conditions 
Op,/Oy = Oat Y = 0 
since 7, = 0 there and pj) > Oasy— o~. Let 
a=x+y/a, B=x — /ao (4.4) 
when 
10°), 0ad0B)] = —- poly } [(a — B) 2]ao}ve"(8) 
so that 
a es "(8 R(y)dp (4.5 
™ 0 Ue (P) (wy )¢ (4.9) 
op EN Sita—p/2ie 
nd 
a Op; 7 ; 
Pi - — ag + Ila) (4.6) 
Js OP 
where f is at present unknown. However, since 
op, /dy = Oat y = 0 
(4.7) 
Op,;/Oa = OP,/OB ata = B 
and hence 
f(a = povr "(a) | R(y) dy (4.8) 
7 0 


The value of /(@) is now found by noting that p, > 0 as 


y— ~ so that on the plate 
Pi = por’ (a) [ R(y) dy (4.9) 
ZW 
Near y rn, R~y and near y = 0, R~y 
so that the integral in Eq. (4.9) is convergent. Further, 
when x is large vg(v) = 0.271 M?2U[v/(Ux)] ?[1 + 
o(1)], so that Eq. (4.9) becomes 
- 7 v l i 
—0.197 po MM P R(y) dy (4.10) 
Ux Ax J 0 


We may now check that all our assumptions are justi- 
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fied when x is large, and we conclude that the fourth 
term of Eq. (3.24) has a contribution from the entropy 
variation in the entire inviscid field. It cannot there 
fore be calculated without knowing the complete solu 


tion. 
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Methods tor Computing Supersonic 
Aerodynamic Flutter Coefficients for 
Two-Dimensional Flow 


ERNEST C. KENNEDY* 


Convair, A Diviston of General Dynamics Corporation, Daingerfield, Texas 


SUMMARY (—), («) = dash over symbol and asterisk denote re 

mm . ee ; ? , imaginary parts of a complex quantity 
The object of this paper is to collect and discuss some of the 
existing methods of calculating supersonic aerodynamic flutter 


coeflicients and to exhibit a new method for evaluating certain INTRODUCTION 


coefficients are expressed in terms of Mach Number and reduced 


associated integrals. The eight real and imaginary parts of the ; , 
, - . ' Se" THE ADVENT of the high-speed plane, some 


frequency (w) by means of convergent series involving the sine intensive work has been done on the theory of 
and cosine of w. It is thought that this paper will be of value to supersonic flutter analysis. Pioneers in this field were 


engineers engaged in supersonic flutter analysis Possio! and Borbely.?. In 1936, the former set up and 
solved the linearized partial differential equation of 
NOMENCLATURE motion of an airfoil oscillating in a supersonic air 


semichord (ft stream. Six years later, Eorbely tackled the same 


problem and arrived at the same differential equation. 


J = standard Bessel function of order x 

M = Mach Number Neither of these two men performed any numerical 
V = free-stream velocity (ft. per sec.) calculations in their papers. Both made simplifying 
) = sin! (1/M assumptions which cause the theory to break down wher 
d = identifying subscript; 0, 1, 2, 3 the Mach Number is near one. 

Q = circular flutter frequency (rad. per sec.) Borbely established the following relationships, taken 


Ww = reduced frequency = 20h/(V cos? @) from references BS and e 


wtf, (M, w) = | ew T)(u/M) du, = £23 
JV 


) ° ° Xr 
= cos uu Jo(uM) du — 1 | sin wu” Jy(u/ AT) du 
7 /J 0 


M? w COS w 
f, (M7, w) = a Ji\(w/M) + (1 — A) cos wly(w/M) + 
w(.M? — 1) I M 
w sin wlp(w/M) — (1 — 2ddofy—1 (WM, w) + (1 — ASP fo MV, ») | + . 
tlw cos wly(w/M) + (w/M) sin whi(w/ A) — (1 — AX) sin wJo(w/ AD) + 
= . ‘ { 
(1 — 2N)wfy i(M, ®) + (1 = o (M, @) ( 
; [w/2-(.M2 — 1)/A?)" Z 
(M,w)=e ™ > A, [Jn(w) + tJ nt+1(w) ], where A, = 3 
n=0 n! (2n + 1) 
and showed that wf) is the fundamental integral used can be expressed in terms of fj, w, and Jf. In this 
in the evaluation of the forces and moments involved. paper no attempt is made to develop the theory or to 
Hence all of the eight aerodynamic lift and moment derive the fundamental equations of flutter analysis. 
coefficients, as defined by Eqs. (49)—(52) of reference 3, The method developed by the writer was used in cal 


culating the table of aerodynamic force and moment 
coefficients listed in reference 5. These eight flutter 
coefficients were computed to high accuracy over a wid 


Received June 12, 1954. Revised and received, October 7, 
1954. 
* Senior Research Engineer. 


310 





Veth 


Ob 


sults 1 
nd th 
f mot 
lapte 


ms. 


Metho 


Bart 
ilua 


le. 


useful 1 
toIEN 

Che 
to that 


Methoc 


isting n 
evaluat 
As i1 


ind illu 


COS 


some 
Dry of 
1 were 
ip and 
ion of 
iC air 
same 
ation. 
lerical 
jifving 


. when 


taken 


1 this 
or to 
SiS. 

n cal 
yment 
lutter 
wide 





range of values of both Mach Number and frequency 


narameter and are available in the above publication. 


\IETHODS OF EVALUATING THE AERODYNAMIC 


COEFFICIENTS 


Vethod of Reference 3 


Obviously, the f-functions of Eqs. (1), 


complex and may be written as /) f, + af,*, 
vhere the dash and the asterisk indicate real and 
ywinary The 


nd /,* can be expressed in terms of w, J/, and the 


components, respectively. functions 


s of lower order by means of the recurrence formula 


[he method of procedure is to first find fy and f,* (by 
breaking Eq. (3) up into its real and imaginary parts 
ind then get /,, (A The 


infinite series (3) converges rapidly for w small and .\/ 


1, 2,3), by applying Eq. (2 
- i 

near unity However for such values of the two param 

factor J/? w?(.\/ | i 


large and consequently these equations give poor re 


eters, the occurring in (2) is 


ults unless each term is evaluated to high accuracy 


s 


nd this is a little awkward to do when w .J/ is a number 
f more than three decimals. This method is not well 
lapted to IBM procedure, partly because of the J 


terms. 


Method of Reference 4 


Barton and Poindexter? devised a neat method for 


evaluating {, which does not necessitate the use of a J/, 
They 
tions and the Fessel functions into infinite series and 


1 


expanded both the trigonometric fune 


takin 
idle. | 


ntegrated termwise. The resulting double series is 
one a finite series in .J/ and 
the other a fixed .1/. 


series converge rapidly for .1/ > 1 and w small, but if 


broken up into two series 
an infinite series in w for These 
is large (>3, say) the convergence may be slow. 
Each /, is obtained independent of fj. The method is 
useful in many cases, but it is not especially well adapted 
toIEM procedure. 

lhe scheme described in reference 6 is very similar 


to that given in reference 4 six years earlier. 


Method of the Author 

In attempting to avoid some of the | mitations of ex- 
isting methods, the writer used the following scheme for 
evaluating the /,'s. 

As in the preceding methods we write 


\ 
W I COS 111 ‘ Ju(u/M) du } ) 


/ 


ind illustrate the method by evaluating f,. Thus 
COs uJ u AJ) du _ 
[ [ 4 VJ (7t AJ)3 (Ut WV ; 2 
— ry = | du 
64 304 


$3402 
6.7816 


= §. 7816 


1.7095 
2.4928 
9 3859 
2 S469 
7.2422 
1. 4963 
2 5978 

8429 
$ GO16 
5 1462 
5.3186 
$ HO09 

5500 


> 
2.4585 


PABLE | 

2 n/2 

, 
were. } 
S402S 6 
S402S 
02797 10) w 
O7549 12 
HHORG 1S 
03301 17 
5S48') 20) 
B4500 22) w 
02772 25 
O3509 IS 
H2O8S8 31 
9IS19 34 
14355 vi 
03424 HH) 
79802 45 
O4018 $6) w 
H5010 14 
S590 53 
21278 jt 
se IS 59 
77205 63) w 
5OSS81 66 
10 


Making use of the reduction formula 


/7 


yi “u- COS U au ul usin u m COS u 
m (m l a u’ 
we get 
Pw 
| cos u du sin @ 
J0 
eas 
] u cos u du w SIN w + COS w l 
« ) 
mas 
rE “u-cos udu 
eo 
W@W) W sin w + 2cos w 2] 


/ u*® cos u du 
if] 


w*lw sin w + 3 cos w 6] 
Pw 
/ u” cos u du 
e/D 
w’ w sin w + # COS w n(n 1) J 


-cos u du 


()) 


Only the even subscripts are used in evaluating fo, but 


the odd ones are included here for later use. 


f D D./ MP? 


where 


D, = I Qo = 


w SIN w 


w SIN w 


+ Dy/ M D¢/ MM 


+ 2 cos w O0.5D 


+- $| cos w 


w Sin w + 2 COs Ww 


D,,-», (n even 


- 0.75Ds. 


Whence 
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where the coefficients a,, = 2" [(n/2)!]—*w"—? are 
listed in Table 1. Note that D,, 17" is not the mth term 


This should cause no confusion. If w 


to 


of series (S). 
is large it is usually convenient to transform Table 1 by 
setting w = 10w’, taking care to shift the decimal point 
properly. 

A little study shows that the series (8) converges for 
all values of w and 7 in which we are interested. It 
converges very rapidly for w small and .J/ large (e.g., 
forw = 1 and J/ = 2, five terms give f) correct to nine 
decimals: for w = 10 and JJ = 2, thirteen terms are 
necessary for the same accuracy). 

Similarly, it can be shown that 


fi = Ey — E3/M? + E,/M* — E,/M®... (10) 


where 


Q,-3(w sin w + n cos w) — 


[n/(n — 1)]E,~», (n odd) (11) 
fo = Bo — B,/M? + Be/M* — Bs/M®... (12) 
where 
B, = d,-4 (w sin w + 2 cos w) — 
= By -2 (n even) (13) 


(n — 2)? 


fs = Gs “ Gs M? oe G7 MW = Go ‘M6 oe (1 }) 


where 
Gn = Qy-5 (w sin w + m COS w) — 
nin — 1) _ as 
= iiy-2 (n odd) (15) 
(n — 3)? 


In using the above formulas note that D,, E,, B,, 
and G, are all functions of w only and that all of the a; co- 
efficients have an even subscript, regardless of whether 
n is even or odd. 

To get the imaginary components—the /*’s—we 


write 
wfo* = f sin uJy(u/M) du = 
0 
f sin u[1 — u?/4M? +u1/64M*...] du (16) 
JV 


Using the reduction formula 


fru” sin udu = u"— [m sin u — u cos uj — 


m(m — 1) fu™-?sin udu (17) 


and proceeding in a manner similar to that described 
above we get 


fo* = Po — Po/M? + Pi/M* — Po/M*... (18) 
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P, = Q,-2 (n sin w — wCcosw) — 
(n — 1) 
Fa (m even 19 
n 


fi® = R, — R3/M? + Rs/M* — R7/M®... (4% 


i = d,-3 (N Sin w — wCcOosw) — 
n 
ae (nm odd 9] 
n— |] 
fo*® = Qe — Qa/M? + QOo/M* — QOs/M*... (29 
Qn = Qy4(" Sin w — w Cos w) — 
n(n — 1) 
(= (m even 22 
(n — 1)? 
fs* = S3 aici Ss M- + S; 1) saad So M* ee 24 
S, = Gy-5 (n Sin w — wCcosw) — 
nin — 1). 
ae (n odd) (25 
(in — 3)* 


Note that if w is very small (<0.1, say), it is neces 
sary to evaluate /,, /2, /3, quite accurately. 
in this case very few terms of the f-series will be needed 
always recalling that \/ 


However 


regardless of the value of .1/ 
is larger than one, say / > 1.1. 


CONCLUSION 


The writer's method is especially useful when only a 
few values of w and many values of 7 are concerned 
It is applicable for either large or small w and each 


is evaluated independently of fp. Since most of the 





calculations for each f, involve a recurrence formula, 
the scheme works quite smoothly on a CPC-IBM | 
machine using an eight-place floating decimal board 
No tables of Bessel functions are required, but accu- 
rate trigonometric tables are necessary. Of course, if 
an IBM machine is available the trigonometric functions 


may be generated as needed. 
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Effect of Rotatory Inertia and Shear on 
Maximum Strain in Cantilever Impact 
Excitation 


MARCEL K. NEWMAN* 
Unwersity of Notre Dame 


ABSTRACT 
A solution of the Timoshenko beam equation is given for 
boundary and initial conditions similar to those obtaining in the 
a cantilevered beam. Curves are drawn for the 
A and for A (w/p, 
These show that for extremely short 


drop test of 
umplification factor as a function of the rela- 
tive impact duration w/p; 
juration disturbances the effects of rotatory inertia and shear are 
of the same order of magnitude as those obtained from the 
Bernoulli-Euler theory with the addition of viscous internal 


damping on the constant Q hypothesis (Sezawa’s equation 


INTRODUCTION 


jibe PROBLEM OF RESPONSE OF BEAMS to non- 
periodic excitation has received considerable atten- 
tion. In 1944, Symonds and Vigness! discussed effects 
of impact duration on the stresses at the root of a can- 
tilever beam. The solution given was based on the 
Bernoulli-Euler equation. 

Theory showed that when an elastic system is sub- 
jected to a nonperiodic disturbance, the response of the 
system is determined by the ratios of the frequency 
fundamental 
If the 
ratio of the dominant frequency component of the 
disturbance to the lowest natural frequency of the 
system is small, vibrations of the system at any of its 


components of the disturbance to the 


frequency of natural vibrations of the beam. 


natural frequencies are not excited appreciably, and the 
response may be found as if the excitation were static. 
As the dominant frequency component of the disturb- 
ance is increased and allowed to become equal to or 
greater than the fundamental frequency /p; of the 
system, the fundamental mode of vibration of the 
system is excited, and this must be taken into account 
in calculating the maximum strains at the root of the 
cantilever. As the pulse duration is decreased still 
further, the Bernoulli-Euler theory predicts maximum 
strain at the clamped end which eventually increases 
without limit. 

The defect in the results of the Eernoulli-Euler theory 
was of course attributed immediately to the omission 
of the effects of rotatory inertia and shear and also to 
the omission of damping. 

In 1947, Mindlin, Stubner, and Cooper” discussed 
results of tests and also results of calculations using 
the Bernoulli-Euler theory, including viscous internal 
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~ 


damping on the constant Q hypothesis, and concluded 


that the 
Bernoulli-Euler equation for extremely short duration 


excessive strain predicted by the simple 
impact excitation was due mainly to neglecting the 
effect of damping which normally provides a mechanism 
for the rapid dissipation of energy, particularly in the 
higher modes. This is in agreement with the sugges 
tions made originally by Symonds and Vigness. ! 

The value of Q used in this case was QO lOO. It 
was found by observations on the decay of free vibra 
tions of the first mode of cantilever specimens. 

It is, therefore, of some interest to obtain solutions of 
the Timoshenko equation,® * which include the effects 
of rotatory inertia and shear with precisely the same or 
equivalent boundary conditions as those offered by 
Symonds and Vigness.! 

There have been given a number of other solutions 
of the Timoshenko equation. Uflvand® treated the 
case of an infinite beam subjected to a concentrated 
load with a function-time history. Dengler, 
Goland, and Wickersham‘ 


using the unit impulse as the time function of dis 


step 
solved the same problem 


turbance. 

In addition, a few studies have been made relating 
specifically to the effect of rotatory inertia and shear 
on the propagation of waves and on the natural fre 
quencies of beams. Davies’ considered the effect of 
rotatory inertia and shear on the wave and group ve- 
locities of flexural vibrations as they compare to those 
of dilatational disturbances. Ayre and Jacobsen* and 
Kruszewski® discussed effects of transverse shear and 
rotatory inertia on the natural frequencies of uniform 
beams. Recently, Schirmer!’ considered the problem 
of wave propagation of flexural vibrations in beams gov- 
erned by the Timoshenko equation by analyzing sep 
arately the shear motion and the rotation. 

More recently, Anderson'* presented a general solu- 
tion of the Timoshenko equation in terms of generalized 
orthogonal coordinates representing bending and shear 
He applied the solution to the case of a 
the midpoint of a 


deflections. 


concentrated transient force at 
hinged beam. 

Also, Miklowitz'* has given a general solution similar 
to that of Anderson. He also made comparisons with 
the Dengler and Goland? solution and with the Uflyand® 
number of 


solution and derived transforms for a 


boundary conditions both for finite and infinite beams. 








0” 


The object of this paper is to present a solution of 
the Timoshenko’ equation for the case of half-period 
sine excitation applied to the root of a cantilever beam. 
This is the type of excitation normally considered typi- 
cal of the drop test. The appropriate boundary condi- 
tions are obtained by means of the use of a variational 
principle. They are equivalent to those used by 
Symonds and Vigness' and Mindlin, Stubner, and 
Cooper.’ 

Results are plotted in a form similar to that used by 
earlier investigators. Amplification factor A for impact 
response is defined as the ratio of the maximum dy- 
nainic strain at the clamped end to the maximum static 
strain at the clamped end which would be obtained if 
the maximum acceleration were reached very slowly, 
as in a centrifuge test. Amplification factor is plotted 
versus frequency ratio p;/w, where /; is the fundamental 
frequency of the bar and w is the dominant frequency 
of the disturbance. 

Another set of coordinates is presented in order to 
make comparisons with earlier results. The ordinate 
A (w/p;) is plotted as function of wp; for two cases, that 
of a long bar and that of a short bar. 

It will be shown that rotatory inertia and shear will 
lead to slightly different limits of maximum strain at 
the clamped end in the case of extremely short duration 
disturbances applied to slender bars as compared to 
thick bars. The distinction is made by means of non- 
dimensional ratio LA, where L is the length of the bar 
and A is the radius of gyration of the section with 
respect to the centroidal axis. All results are given in 
nondimensional form when such form appears to be 
indicated. 

Specifically, numerical computations were carried 
out for the thick bar with L/A = 30 and for the slender 
bar with L/A = 300. These values were chosen with 
the intention to cover the approximate range of prac- 
tical significance of values of L/K. 

It was found that for short duration impact, the value 
of the ordinate A (w/p;) approaches a limit of 4.89 in the 
case of the slender bar and a limit of 4.08 in the case 
of the thick bar. These values compare to that re- 
ported by Mindlin, Stubner, and Cooper? in the case 
of a bar computed on the Bernoulli-Euler theory but 
with internal viscous damping according to the constant 
Q hypothesis. This value reported is A(w/p;) = 3.55. 
It shows that effects of rotatory inertia and shear are 
of the same order of magnitude as those of constant Q 


damping. 


(1) EQUATIONS OF MOTION AND BOUNDARY CONDITIONS 
FOR CANTILEVER BAR 


The equations of motion for the elastic system sub- 
jected to a drop test and the appropriate boundary 
conditions arising from consideration of rotatory in- 
ertia and shear are obtained conveniently by minimizing 


the integral 


ff’ (T — U)dt (1) 
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where 7’ — U is the Lagrangian function, 7° is th, 
kinetic energy, and l’ is the potential energy of th 
beam. 

The generalized coordinates chosen are w, the trans 
verse displacement of the beam element, and ¥, th 
angle of rotation of the element due to bending. Th, 
slope of the neutral axis of the beam is 

w/dox = V+ 4 


where x is the coordinate parallel to the axis of the bea 
with origin at the free end and y is the shearing strain 
The bending moment is 


M = —EI(OW/oOx) 3 


where /:/ is the flexural rigidity of the beam and OW d 


is a first approximation to the curvature for small 


values of WV. 
The shearing force is 


Q = K’AG y = K’AG|(Ow/dx) — V] 


where A’ is a numerical constant depending on th 


{ 


shape of the cross section and on Poisson's ratio," 4 
is the cross-sectional area, and G is the shear modulus. 


The kinetic energy is 


ae *! /ow\? 
i= > pA f (**) dx + 


bas * (aw? 
K? pA } ( ) dx (5 
2 J0 Of 
The potential energy is 


1. f* Sow’? 
= KI dx + 
2 J0 Ox 
] ac Ww 
K’AG } ( — ¥) dx (6 
2 J0 Ox 


where A = radius of gyration with respect to horizontal 


principal axis through centroid of cross section of bean 
and p = mass density of beam material. 

Definitions, Eqs. (5) and (6), exhibit the degree o 
approximation intended. 
indicated, Euler's differential equations of motion ar 


obtained, together with the natural boundary condi- 


tions that must be satisfied in order to minimize in 


tegral (1)—namely, 


— pA (0°w Ot?) + 
K’AG|(0?w/Ox?) — (OW/Ox)] = 0 (i 


—K*pA(O°W/Ol?) + EI(O0?W Ox?) + 
K’AG[(Ow/dx) — VW] = 0 


Eqs. (7) and (8) can be combined by eliminating th 
generalized coordinate V. The result is the well-know! 
equation 


O*w p(l +e) OW 


ox! EK dxdt? | 
pe OfwW 4 p wW_o 
E? of K°E Of? 
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The natural conditions at the ends x = O and x = L 
ire 

(OV Ox )db¥ = 0 (10 

[(Ow/Ox) — V]éw = O (11) 


where L is the length of the beam. 
For arbitrary 6W the condition of zero moment is ob- 


tained 
dy /dox = 0 (12) 
Substitution in Eq. (7) leads to the form 
O°w/ Ox?) — (pe/E) (0°w/dl?) = O (13) 
If 6w is arbitrary, the condition of zero shear is ob- 
tained 
[((Ow/Ox) — VW] = O (15) 
Substituting this in Eq. (8) and using Eq. (7), the 
condition of zero shear takes the form 
O*%w Ox*) — [p(1 + €)/F] (0%w/dxdt?) = 0 (16) 
If the bending moment is arbitrary, Eq. (10) requires 
that 
ov = 0 (17) 
and for arbitrary shear, Eq. (11) requires that 
6w = 0 (18) 


The quantity ¢€ is a dimensionless constant defined as 


follows 
e= E/(K’'G) = 2(1 + »)/K’ 


where » is Poisson’s ratio. 

eis adjusted in accordance with the suggestion made 
by Mindlin'* so that the limiting transmission velocity 
for extremely short flexural waves approaches that of 
Rayleigh surface waves. 


For a wave solution 


Eq. (9) takes the form 
ct — (c,7/K2y?)c? — (c,2 + c,2)c? + c/?c,2 = 0 (14 


where c, is the velocity of long extensional waves 
= E/p; c.2? = K'G/p; c;7/c, = E/(K'G) = e«: 
u = 27 X; and J is the wave length. 
Eq. (14) implies two modes of transmission’" 
given by 


a | l 

2% (1 T Ky? + ? 2e| x 
1 4 ; /( ean Et 
yi 2 \ e/\ + Kz? §) \f (14a) 


Confining attention to the lower transmission mode, 
the short wave-length limit with 1 (K2u?) = O will be 
given by 


c? = K’'G/p 


T 


I 


L 
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This corresponds to Mindlin’s'* equation c? = A°(G ‘p). 
For materials with Poisson's ratio vy = 0.2, A? will be 
0.828 and e will be 2.9. For v = 0.3, «€ = 3 


(2) CANTILEVER BAR IN Drop-Test EXCITATION 

The initial conditions and end conditions for the case 
of excitation in a drop test of a cantilever bar of length 
Land with origin of coordinates at the free end are 


[zw], 90 = 0 (19a 


[Ow Wt),-0 = —V (19b) 


v (1 = 3 u (: — "| (19¢ 
\ Ww Ww 


[v].,-, = 0 (19d 
After some manipulation, the condition of zero angle 


of rotation V at the clamped end of the bar x = L takes 
the form 


( 37) € Ow ] yra ) 
oo a oe 0 (19e 
Ox? Cc,” Oxo? eK* Ox f eas 


The boundary conditions at the free end of the bar are 


( “Ww € =) 7 
— — = 0 (194 
Ox? cr 7. =6 


| 37y) ( l + €) S7y | 
ox* ¢; e@xot i. 


In the foregoing equations u is the unit step function 


| | 


0 (19¢ 


u[(r/w) — t] = 1 for [(r/w) — t] > 0 
u[(r/w) — t] = O for [(r/w) — t] < 0 


| is the initial velocity of the cantilever at the instant 
of impact. J is of such magnitude that the maximum 
acceleration I’w is large in comparison with the acceler 
ation of gravity so that the latter may be neglected. 
This initial velocity corresponds to an initial kinetic 
energy storage in the beam of !/. m\*, where m is the 
mass of the beam. 

Solution of this problem is accomplished readily by 
use of the Laplace transformation 

Y = { e ” w(x, tdt 20 
wv V 
where p is a complex parameter. 

The partial differential [Eq. (9)] and the conditions 
[Eqs. (19a) to (19g) ] are thus converted to the following 
equationsin }.. 
aey 1 + « m {2 ¢ 


— b= = os 


dx* c I" dx . 


ep* . b= l 
(:? + P )r--(% 4 - )r (21 
c;* K°c,;? Cc; K°c,? 
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Atx = 0 (1 + 6) (l+6). As 
: as Gs FF th 
(d?Y/dx?) — [(€/c,?)p? VY] = (€/c,7) I (22a) a ae 
$17 /3.28 ei Si ; eae | te 1c," 
(d®Y/dx*) — [1 + ©)/c?]p?(dY/dx) = 0 (22b) ¥ = = dias V(a 
(1 + e)? (1 + €)*A2p sa lud 
! ude 
Atx =L a 
It will be noted that y,° and yu? are of opposite sig -) 
_— 7o—(pPu/w) §,, 2(p2 ia re 8 
7 orm 1@/ [p?(p? + w*) |} rhis means that the corresponding values of ,; and , It 
[V/(p? + w?)] (22c) are at quadrature in the complex plane. Becaus¢ wi 
ree a - ; > symmetry of Eq. (23), one of the ble signs j ton 
(d3¥ /dx3) — [(ep?/c2) — (1/eK2)] (dY/dx) = 0 (22d) the nmetry o q. (23), one o the double sign we 
Eqs. (24a) and (24b) is redundant. They are dropp 
The complementary function of Eq. (21) can be since they will contribute no new solutions in Eq. (2 Eq 
written in the form Eqs. (24a) and (24b) are therefore given as follows ue 
Y = ¢ sinh wx + c@ cosh wx + 1+ .e l+e ns 
. -= a h2 x W(X, | 
C3 SIM pox + Cy COS pox (23) ae 2c)? f 2c; t 
provided the following quartic equations in 4, w2, and le ke 
p are satisfied. (1 + ©)? (1 + 6)*K%: 
' ee rh 
mit — [C1 + ©€)/c;7] pu? + (epic) + Ite. oe ; 
[p? (K°c,7)] = 0 (28a) ae ee ae ee de 
tL 4 °C; = 20a 2c;? 2c/? Acubl 
and be 1c,” exist 
— —— 24d 
ore or 2 (1 + 6&)*K’p? ~ An 
pot + [(1 + ©€)/c;7]p7u2? + 
(epi/c,t + [p?/(K2c?)] = 0 (23b) The signs of w,* and p.? in Eqs. (24c) and (24d) r 
flect the fact that Eq. (23) is chosen hyperbolic in u shows 
In general, each quartic has two roots in the squares and harmonic in w2x for positive real values of 4; and y Thev 
of the parameters y; and yo, respectively. These roots The solution of the differential Eq. (21) with bound lee 
are ary conditions (22a) to (22d) leads to the Lapla corres 
transform }(x, p) of the solution w(x, /). » gov 
ian tts, G49. / go\ 
2 — 2 + 2 ete a : } 
My 2e;* p 2¢,2 p VoX(1 + e~°*/*) (x, p) | that ¢ 
V(x, p) = (p2 + w)ACL oe ae tion 2 
l a de," 4 ee nee , Eq 
~ - —— (24a) 4. 
(1 + e)? (1 + ©)?K2p? where axis. 
Int: 
P(x, Pp) = mime! (c?/ AK?) [cosh w(L — x) + cos po (L —x))] + [(e — 1) (ui2c? — ep?) + (c?/K?)] cosh mL cos pox + 
[—(e — 1) (ue2c,? + ep?) + (c,2/K?)] cos woL cosh mx} - 
(uo*c,*/AK?) sinh wZ sin pox + (ui7c,2/ A?) sin wel sinh wx (2 and 
and 
A(L, p) = mime} (2c;2/K?) + [(2¢,/K*) — (<e — 1)2p?] cosh mL cos wel} + [(1 + ©) p?/ A?) sinh wZ sin wl (27 ie 
| the ec 
venier 
| 
It can be shown that wo = 1Mp’[p + N(p? — a?*)” 24i) | 
‘ 
#(L, p) = ACL, p) where AL, 
The use of two parameters yw; and uw, in the Laplace M = [(1 + €)/2c,?] 
transform V(x, p) [Eq. (25)] insures that V(x, p) ; 
a re ° ° : ° N = (1 — 6)/(1 + @) 
[Eq. (25)] is a single valued function of p at the branch 
points of wu; and po. a = 2¢,/[(1 — eK] | where 
The branch points of uw; are a ‘ ‘ : ae - 
- he point p = +2(c,/A Ve), which is a zero of tl 
p= 0, p = +2c¢,/[(1 — 6)K] function 


The branch points of ye are z= [p + N(p? — a?)] 


p = 0, p = +2c¢,/[(1 — 6€)K] is not a branch point of 4; and we. This is shown by pel 
fice a ‘ forming a circuit in the complex plane around this 
This can be shown by writing yu; and yw in the form : : é ‘ae te 

point and noting upon completing the circuit that 5! 


mw = Mp [p — N(p? — a?) ] (24e) single valued. 
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A single circuit in the complex plane around any one 


of the branch points of uw; and yw either changes sign or 
changes the value of wu and yw» in such a way that 
») [Eq. (25)] remains single valued. It is con- 


= J 


cluded that the only singularities of ](x, p) [Eq. 


95) |are poles. 
It remains to find the solution w(x, ¢) of the problem, 
p) is the Laplace transform. This solu- 


of which } (x, 
Since the in- 


tion is the Laplace Inverse of ](x, p). 
version process is unique, w(x, f) will be the solution of 
Eq. (9). 


a 


the contour integral 


ul 


rhe inversion is accomplished by means of 


ly ¢ el +e ‘°"* )b(x,pidp 
aici ot B p?(p? + w*)A(L, p) 


vg “7 (28) 
271 Br, p* 


The solution is the sum of the residues at the poles 
f{ the integrands of Eq. (28). In addition to the 
double pole at = O and to the poles p = +7, there 
exist zeros of A(L, p). 
An analysis of the equation 
A(L, mL, po, p) = 0 (29) 
shows that these roots are uw = yw), and wl = po,L. 
They represent the characteristic values of the prob- 
The values yu, and yw, occur in a one-to-one 
They are related to a set of values of 
The values of p 


lem. 
correspondence. 
p governed by Eqs. (24c) and (24d). 
that correspond to the roots u,/ and po,L of the equa- 
tion A(L, wil, wel, p) = 0 are poles of the integrand of 
Eq. (28). These poles are located on the imaginary 
axis. They represent the natural frequencies p = 
+ip, of the modes of free vibrations of the bar. 
Introducing now the circular frequency 

w = C/a (Ky €) 
and 

Pn = C10n/(KY €) 


where a and a@, are real positive numbers 0 < a < @, 
may be written in a con- 


the equations for w,2 and p2? 
venient nondimensional form. 


A, (Mins Mons Pn) ((2€7 
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If p = w = 1t|(c,a)/(Kve)] is substituted in Eqs. 
(24c) and (24d) the following is obtained 
(l + « (l + «€ 
K Ma” = = - a? -— . x 
aw€ ow € 
la = £€ le 
a T 242 
(1 + «€ 1 + €)"a 
(1 + e« ; (1 + e« 
K M = = x 
2e Ye 


ie le 
a- | T (24h 
Vi + €)? (i + €)*a" 


Ky, will be real for 0 < @ < 1 and imaginary for | 
a<o,. Ky willbe real for all values of a, 0 < 

The by computation of 
A(L, ml, wl, ta), using appropriate values of wL, 
All roots in p of A(L, wsL, pol, ta 0 are 


< x 


SS = 


roots were obtained 
pol, and Qa. 
simple poles. 

The first six roots are listed in Table | as function of 
a, for the short beam, L A = 30, and for the long 
beam, L/A = 300. 

The solution is written for the period of the dis- 
turbance 0 < ¢ < w/w and for the time after the dis- 
turbance has ceased to exist, rw < ft < «a 

For the period 0 < ¢t < ma, the displacement is 


| i P(x, 1w } 


Ww Als. lw) 


sin wi + 2Vw? X 


Wx,Wy=- 


P(x, ip,) sin Ppl 


a b,*(p,* — w)A’(L, 1Pn) 


(30) 


For rw <t< @ 


Tv , 
(1 2) = are x 
Ww 


P(x, tp,) cos [(rp,)/(2w)] sin p,[t — (@ 


n Prn?(Pn? — w?)A'(L, tpn) 


W(X, ti) = 


2w) | 


1 
(91 


where w = frequency of impact and p, = natural fre- 
quency of mth mode of vibration. A(L, iw) is obtained 
from Eq. (27) by substituting p = 1; (x, tw) is ob- 
tained from Eq. (26) by substituting p = iw; ®(x, 7p,.) 
is obtained from Eq. (26) by substituting p = ip,; and 
A’(L, tp,) is the derivative of A(L, p) with respect to 


p for values of p = ip,. It takes the form 


K*) + [(e — 1)*p,? + (2c,;?, K?)] cosh w,Z cos w2,L} — 


2ui,m2,(€ — 1)°p, cosh w,L cos po, + [2(1 + ©€)/K?]p, sinh w,Z sin we, + 


{ Minkton[(e — 1)2p,? + (2c,?/K*)]B, — 
{u41,M@2,[(€ — 1)*p,? + 


where 1), is the value of u; for the nth mode, po, 


[(1 + €)p,p?/K7|C,} sinh w,Z cos wo,L — 


(2c;* K?)]C,, + [1 + €)p,” K?]B,} cosh y;,Z sin po, L 


is 27/\,, proportional to the wave number of the mth mode. 


2 , on 2. an u1,Pn 
A,| , Me ee: ae a 4. be, f : 
(ins B2nr Pn Pr Ky, [2u1,,2¢,27 + (1 + €)p,?7] Kus, [2ue,2c:2 — (1 + ©)p,?] 
. - n 
B, = limL— =L . Sy = - . ~ q 
p=pn Op | Pn K au (2:2; + + e)pn?]f 
he i “ J 2n Pn 
C, = lim L “se | 
; — b | Pn ¥ Tin Mate? — 1 4 p14 
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Solutions (50) and (31) satisfy the differential Eq. nt : 
ee ait ee _— yhe 
(9) and the initial conditions and boundary conditions Characteristic Values for Cantilevey 
Q: (19¢@) inelusive : 
(19a) to (19g) inclusive. L/K = 30 short beam L/K = 300 long bea 
a Mi iB Mon a, ub / im 
(3) THe AMPLIFICATION FACTOR 0 0064242 1.8336 1.881 0.00006720 1.875 1 .& P 
0 .0390157 $ 446 t 644 0. 00042708 } 6X8 1 RoC 
The results obtained may be expressed conveniently 0.1014267 6.864 7.743 0.00116583 7.860 7g 
‘ E ‘. pied sag a 3 ; 0.1819940 8.760 10.791 000224506 10.860 10 4p: O24 
in the form of an amplification factor for strain at the 0.273441 10.029 13.836 0.00375991 13.950 14 146 
=9947 nie baer seo a it 
clamped end of the cantilever. It is defined as the ratio 0.323474 10.890 16.875 0.005950713 17.010 17,99 
of the extreme dynamic fiber strain at the clamped end 
of the bar to the extreme fiber strain at the same loca- acteristic dynamic behavior in the range of short dy 2 
tion that would result if the maximum acceleration ration impacts, it is preferred to discuss the quantit 
lw were reached very slowly, as in a centrifuge test. A(w p;) as function of wp). This quantity is propor 
The extreme fiber strain at the clamped end of the tional to the maximum dynamic strain at the clampe 
bar is given by end 
&. = d/ 20 Ox) 1 (32) A(w pi) =  ( % e/a) (K2/L? (le | 
where d is the depth of the beam. where 
The extreme fiber strain, at the clamped end, that nies . 
ae , : i a? = EK? p = K* whet 
would result if the maximum acceleration Jw were 
reached very slowly, as in a centrifuge test, is a, is defined by fundamental frequency of bar as 
€,, = (VwLh?d)/(4c,? K?) (33) pi = (€:a1) (AV) 
The amplification factor for strain at the clamped I” = initial velocity of bar at moment of impact 
end is here defined as the maximum value of €,,/ €,,. d = depth of the bar; and K L nondimensional 
enn ~ — ; atio of radius of gyration of cross sectio1 noth of 
A = §|[2K%c;2(0W/Ox),=1]/(Vo@L?)} maz (34) rat 5. . 10n to length 
bar. 
4 ‘ shoati Pap Py 2 1c > aga: ~ > 
he amplification factor, Eq. (34), is plotted against It will be noted that A (w, p;) depends on (A L 1. 
: me ¥ = se At "= 2); "= 2 : ° o = 1 . 
p:/w in Fig. 1 for values of L/K 0 and L/ A 300. its value will be different for a slender bar than for 
The amplification factor for the Bernoulli-Euler sys- thick bar. In the case of L/K = 39 
tem is also plotted on the same figure. 
or , " ; A(w/p,) = 1.18|(ae,.)/(Vd)]} 26a { 
Differences between the curves for L/AK = 30 and Pi | : sia 
L/K = 300 were found almost imperceptible and and in the case of L/K = 300 of i 
arise primarily in the range of extremely low frequency 
A (w/pi) = 1.13 [(de»)/ (Vd) ] 36b Irequ 


ratios p/w. 
It may be seen that the Bernoulli-Euler theory is A(w/p;) was computed from values of w(x, ), E 
; ; dad dae reaso! 





close to the more exact theory over most of the range of (30), for the period of duration of impact 0 <¢ <7 “Hl 
: , ‘ ae of the 
frequency ratios. However, in the range of extremely and from the values of w(x, ¢), Eq. (31), for the period Dt 
small frequency ratios p,/w, the Bernoulli-Euler theory r/w <i < ©, after the disturbance has ceased. 
predicts infinite strains, not in agreement with the Curves of A(w/p;) were plotted as function of « ; 
more exact theory which takes into account effects of for L/K = 30 and for L/K = 300, a range that covers 
rotatory inertia and shear. This lack of agreement for most practical cases 
extremely short duration impact is not apparent from , ; : ‘ age 
ic. | ‘ Equations for A(w/p,) have the following form. 
Fig. 1. 
This is an inherent shortcoming of the amplification- eee , where 
: pally During impact, forO0 <¢ < 
factor curves. In order to bring out clearly the char- w 
Aft 
w Ka j [0° (ta) /Ox?], 7 -_ 2a[0°@,(ta;)/Ox*,-7 pb 
A = 2Ve— €° _ lf sil wt — = ere sin pit) — 
Pr L?ay | Ai(ta) (a," — a*)A;’ (2a;) 
ee 2 (0°, (za, ) Ox? ], I , _ 
‘s+ <—, » - sin p,t (vi 
J? Qa, 2 An~ ee Aga 
=- = 2) Ay’ Ge,) sa 
a~ 
After impact for r/w <t < @ - 
his v; 
w ~ K? 1 GQ [0°*,(ia,)/Ox?],- 7, cos [(4/2) (a,/a)] . 7 
i ai — — FZ, ee sin p, (¢ — ” 
pi Le wm 1 [1 — (a,?/a?) ]Ai’(ia,) 20 





Mm pact 


sional 


eth o 


1 for 





MAXIMUM STRAIN IN CANTILEVER IMPACT EXCITATION 10 


vhere 
A, (ia (Kyu) (Kye) }2e + [(e 1 )2a? + Pe] cosh wL cos pol} 1 + e)a’? sinh euZ sin wl 
0°@, (fa) ON" },=1 Ky Nye) (1 r € a? + (3e — 1)] cosh wL cos wel} 
; l é €/a + ?!' snh w LZ sin pw L 
02? (7a ad Ku (Ayu (1 + €) F e— l1)-a s€ 1)] cosh w,,Z cos po,L! 
1 + €*) ela 2! sinh w,,Z sin po,L 
) a ( 
A Y Vv i a E = Bi = r 2 | cosh vi L cos v1 Ly 
€ — 
2(e — 1)?a,2 V1 — a,? (1 4 . ae 
cosh y,,Z cos wo,L 4 sinh y,,L sin po, L 
€ V «é 
L j : of a,” ] 3 1 + e€)a (Au : ' 
© Ku V i a (e — 1 + 2] Fr. Cr sinh a L COS UU B 
AK 4 L € ] Ve 
L j a | 1 + e)a Au | 
: (Aw Vv i w* ih te 1) t 21G, 4 F cosh yw, sin po,L 
A { s € | Vv € { 
where 
( Ze 
\a } € + (te a —(Il1_+e)WN l— € + le a } 
) 
aw ) 
if = 
a rt] € T le a |} 7 € \ ] é 4 le a : 
. Ye 
y, 
a fd —_ 2 + (tea 1+ e«) Vv l € + (4e a ; 
. D¢ 
G, 
2.73 [( l —eé« + (dea t l+e) V1 € le a : 


A(w pi) is a continuous function of wp; for the period of impact O tf < wm wand for the period after impact 


rwo<i< oo. Atthe time? 

When w >, Or a a, ~1.e., when the frequency of the half-period of impact is equal to any one of the natural 
Eqs. (37) and (38) cannot be used. The form of these equations is such that cases of ‘‘reso 
There are obvious physical 


= mr w, both Eqs. (37) and (38) assume equal values. 


frequencies of the bar 
nance’ lead to indeterminacies that can be determined by means of a limit calculation. 
reasons why strains in such cases must be finite and determinate. When a, represents any one of the natural modes 
of the bar, the form of Eq. (37) changes as follows. 


During impact, for0 </ < m/w 


( i? a [0° (Za,) Ox* ] | A > [O° 1a Ox" | ] : 
x au 1 cos wh Ven, & | TT Raa 
p aN TL? a@ A,’ (ta, Qa LD? a,*a |(a@ a ] lA, (1a@, 
kK a 2 
Ye sin wf 9) 
L? a 
where 
v (C; Q (Ky t 


After impact, for m/w <i < 


W Kk? | [0° (ta,) Ox? |, / 
A 2Ve a Tt COS wl + 
p LE? a, A) (ta,) 
K? |] > (O°, (7a, ) Ox? |, 7. COS [(a4 2)(a,/ a } Z 1 \ 
SVE, a : ae sin P, \t - (40) 
L’? QA) a,*a, (1 — (a@,°/@ ) JA, (1a, 2w 
Again it easily can be shown that Eqs. (39) and (40) give the same value for A(w p;) att = mw. For 
W@W = (C; Ay) (K V €) 
this value is 
{ Ww a K? (O°? (7a,) Ox? |, i] ; K? l = (0°, (za, Ox? j, 1 - Tp, 1] 
= ey ¢ ° re _ Ve : a : “CF sin 
Pr L? ay A; (1a, L? Q, a,#~a [(a@,,° ot ee 1 JA, (la, w 
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R/W 


Fic. 1. Amplification factor ‘A’’ as a function of p;/w 


The series in Eqs. (37) to (41) appear highly con- 
vergent for small values of w/p;. Convergence is 
satisfactory for extremely large values of w/ py. 

The nature of the convergence of the strains is best 
illustrated in terms of contributions to A(w/p,) by the 
various individual modes. This is shown in Table 2 for 
three typical values of w/p, and for L A = 30. 

It will be noted that one value of w, p; chosen is small 
and another is extremely large. The small value repre- 
sents a condition arising during the half-period of im- 
pact. The first figure listed in that case as contribu- 
tion by the first mode contains the effects of the forced 
response as indicated by Eq. (37). It is the maximum 
value arising due to beats between the forced response 
and the fundamental. 

The final value of A (w//;) in each case is obtained by 
summation of the contributions by the individual 


modes in their proper phase. 


TABLE 2 
Values of A(w/p,) as Contributed by the First Six Modes for 
Three Typical Values of w/p,, L/K = 30 


Mode w/p, = 0.221 w/pi = 10 w/p, = 10° 
1 —(). 25121086 —1.8217899 —1.8260538 

2 +0 00065495 —().9050575 ~(0 .9875643 

3 +(). 00005248 —(). 2829975 —().5853914 

4 +().00001042 —(). 1264207 —0. 3424051 

5 +0  QOO003809 +0 .0123078 —(Q . 2290035 

6 +(). 00000112 +0 .0044955 —0.1543817 
A(w Pp) —(). 25048880 —3.1194623 —4.0747998 


Curves of A(w/p,) were plotted as functions of wp, 
for L/K = 30 and for L/K = 300. In each case the 
first six terms were included in the series, and it was 
found that little was gained by carrying more terms. 
The results are shown in Fig. 2. 

In the case of the short and thick bar L/AK = 30, it 
will be noted that the limiting value of A (w/p,) for large 
values of w/p; is 4.08. Similarly, the limit of A (w/p;) 
for large values of w/p,; in the case of the slender and 
long bar L ‘AK = 300 is 4.89. 

A (w/;) is proportional to the maximum strain e,, at 
the clamped end. In the case of the short bar, L/AK = 
30 
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Ww ty € K* Em ae, 
A = — Se Te = 12a 
pi a [?/] Vd Vd ’ 


In the case of the long bar, L/K = 300 
Ww ty € K? en de, 
A = : a = ie 12h 
pi ay By Vd Vd 
where 
: EK? ae 
a? = = K¢/ 
p 
|” = initial velocity on impact 
d = depth of bar 


From Eq. (42a) the maximum strain for impact of 


extremely short duration at the clamped end of a short 
bar will be 


€m = (4.08/1.18) (Vd/a) = 3.46(Vd/a) (48a 
For the long slender bar, the maximum strain will be 
€m = (4.89/1.13) (Vd/a) = 4.34(Vd/a) (43b 


In the case of the constant Q hypothesis in the 
Bernoulli-Euler theory 


é, = (3.55/1,.188) (Vd/a) = 3.12(Vd/a) (4% 


This shows then that the effect of rotatory inertia 
and shear on the maximum strain at the clamped end 
of a bar in short duration impact is of approximately 
the same order of magnitude as viscous damping 
would be on the constant Q hypothesis. It is instruc- 
tive to note that a thin slender bar will be subjected to 
approximately 25.5 per cent higher strain than the 
thick short bar at the clamped end in short duration im- 


pact. 
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The Buckling of Isosceles Triangular Plates 


HUGH L. COX* ann BERTRAM KLEIN* 
Douglas Aircraft Company, Inc. 


SUMMARY 


Approximate solutions are presented for the elastic critical 
buckling loads of isosceles triangular plates that have either 
lamped or simply supported edges. Loading conditions of 
uniform axial compression along all edges of the plate and of com- 
pression along the base edge and shear along the other two edges 
we treated. The results of the calculations are presented in 
graphical form so that an engineer may determine rapidly the 


lowest critical buckling load 


(1) INTRODUCTION 


Fee ENTLY IN THE ANALYSIS of aircraft structures 
it is necessary to determine the buckling load of a 
triangular plate. Usually the plate is either uniformly 
axially compressed along all edges, Fig. l(a), or com- 
pressed along the base edge with the load taken out in 
shear along the other two edges, Fig. 1(b). The ap- 
proximate solution presented in this paper is based 
essentially upon the method of collocation—i.e., for an 
assumed deflection function, which satisfies the bound- 
ary conditions, the governing differential equation for 
Com- 


parisons are made for certain cases with other known 


the plate is satisfied at a finite number of points. 


solutions, and the agreement appears to be satisfactory 
for most engineering purposes. 
(2a) SIMPLY SUPPORTED PLATES UNDER THRUST 
Let x and y be coordinates in the middle surface of 
the plate as shown in Fig. l(a). The differential equa- 
tion for a uniform plate subjected to thust is 


DA*?A*w = N,(0?w/dx?) + N,(0°w/dy?) (1) 


where 
D = plate stiffness 
w = deflection of plate, positive downward 
N, = load per unit length acting parallel to xv axis 
N, = load per unit length acting parallel to y axis 


The boundary conditions are 


(2 4 9 0 : 
v = (3) 
On? aa oe a ee Ne 


where is the normal direction to the lines 


x = +(a/h)y 


and ¢ is the tangential direction of any line. Along a 


rectilinear edge the term 0°w/Ot? is equal to zero. 


Received July 20, 1954. 


Strength Engineer, Santa Monica Division. 


The deflection function 


w= (a,Fy + asks + aks 


where 
F, = 2 sin? (ry/2h) sin (ry h) 
2 = 2 sin? (ry/2h) sin (2ry/h 
Fz; = 2 sin? (ry/2h) sin (3ary/h 
a = generalized coefficient 
k = odd integer 


satisfies Eq. (2) exactly and the condition 
(O°t0/ Ot") «m= tathy « © 


from Eq. (3) at some point in the region of (h 2) < 


y < (3h/4): the condition that 
O*w 
9 ” 0 
Ov" y=h 
is satisfied exactly only when x = 0. Selected points 


for collocation will be on the line x 0; therefore, no 
serious error should result from not satisfying com 
pletely the boundary conditions. The first term of the 
deflection function given by Eq. (4) appears physically 
to have the shape that one might expect for a triangular 
plate buckled into one wave. 

It will be necessary to differentiate the deflection 
function given in Eq. (4) and substitute the proper de 
rivatives into Eq. (1). 

Let 


F(y) = oF, + ao Fs + @ F; (5 
Eq. (4) may be written as 


ow" = 


w= F(y) cos k(r/2) (h/a) (x/y (6) 


Derivatives of watx = Oare 


O*w wr h\? a 
- = —F(y) (k ) 
Ox" J cmt 2 ay 
O*w i T h\4 
: = F (y) k Ss 
Ons cue ‘ 2 ay 


9 = O°F(y)/Oy? (9) 


(O°w Oy" ly = 


(Oo 40y) oy" lzuQ = 
( ow) Br ( h ~~ 
Ox*0y"/ r=0 : fy" a Oy" 
D2 a2 (*)’ OF(y) 3h? (=)’ ; 
_ ae Fly 
y? a Ov 2y* \a 7 


O'F(y) /Ov' (10 


(11) 
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The derivatives of /(y) are written as 
Fy’ 
Fy" 


If Eqs. (12) are substituted into Eqs. | 
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). | ] 
Til cos TV ia — (06 ica } 
ah)? [—sin (rv h 2 sin (2rv h)] 
(7 h)* [sin (ry h) — S sin (2ryv h)] 
7 2m} ) 1) I Ty 
2 Cos - cos cos 
as h 2 h 2 h 4 
a \* , 2Ry 9 . dry : .. ey 
—-+ sin a sin + we 1 
h h 2 H 2 h | 
‘ ‘ *) a | 1») 
er \" , Sey Sl | sry | . = te 
16 sin —_ sin — 5 sin 
h 2 H 2 1 | 
7 ST) tary 2ry 
3 cos —? COs — Cos 
h 1 1 1 
a _ omy : ! . 2ry 
-9 sin + S sin + 2 sin 
h ay 
a \* , omy . 4try oo ee 
Sl sin — 128 sin — Ssin 
h | i 1 1 


7) through (11) and the resulting derivatives are substituted into Eq. (1), there 


results 
Aa + Ba» + Caz = 0, at y = h/2 
Ja, + Eas + Fa; = 0, at y = 2h/3 13 
Ga, + Ila. + la; = 0, at vy = 3h/4 

where 


h\4 kh \ 1 Ne | 
r(k + mr? + 22.00490 (k — Bik + 1} 
a a t a a 


394.7860 — 21.5740[Rk(A/a) |? — 48 


h\3 _ = h\? h\? ] 
C = —r-(k — 799.4416 — 179.9193 [Rk + 8 k + 9 | 
a a a | 


h\3 
B= —2.34211 (i ) = 
a 


— 592.1790 + 21.784 


- n\* = ie Ne 
41342 + 5.26442 (R —- B| 0.42187 (k + 1.0 
a a 


a? x h\? - 
+ B| 0.42187 (Rk + 2.25 
a a 


S1.42461 — 41.1351S4 (: 


IS[R(A a) |? + 38 


: h\4 7 ay h\? : h\? 
G = 235332 | k + 8$5.93610 + 1.77300 {| R — 3] 0.53649 [R + 2.70711 
a a a 


(ei) 
32809 | Rk — 
a 


and where 


B = Nh?/D 


Values of 3 for various ratios of ha may be determined 


from the condition 


— h\? eo h\? 
144.04924 — 60.45476 [eR + 8B] 0.75871 k + 7.53553 
a a 
na \* 7 hk \? 
— 8} 0.53649 [| kR + 8.36396 
a a 


(15) , B 


~ 
N 

my ee 
II 





o 


here 


BUCKLING OF ISOSCE 


vs in graphical form, on semilog paper, the 


Rio. 2 shov 


between 8 and / a for a practical range 


relationsh}] 


values. For an equilateral triangle, from Eq. 


6), 8 .6. The exact answer! is 6 = 47 The 
hove solution 1s somewhat conservative for an equilat 
eral triangular plate; for the case of / a 1.0, the 


jue of 6 given in Fig. 2 differs by less than one-half 
fone per cent from a value obtained by Wittrick.* In 
veneral, the curve in Fig. 2 affords a rapid approximate 
miswer, Which appears to be on the safe side for plates 


having various ratios of hh a. 


SUPPORTED PLATES UNDER SHEAR AND 


THRUST 


2b) SIMPLY 


For the load conditions shown in Fig. 1(b 


If the first terms in brackets in Eqs. (14) are multi- 


plied by a h)*, the resulting equations will be suffi 
cient for determining the lowest critical value of .V. 
ha. 


.O, the value of 8 given in Fig. 3 differs 


Fig. 3 gives values of @ for different ratios of 


When / ‘a l 
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FIG. 1! LOAD CONDITIONS CONSIDERED FOR 


IANGULAR PLATES 
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by less than 1 per cent from the value obtained by 
Wittrick.% 


Consider the plate shown in Fiz. I(a r 


) 


3a) CLAMPED PLATES UNDER THRUST 


1e be und 


ary conditions are 


IS 
Ow OY) y=4 Ow O71 0 
The deflection function 
rh ay Fy + asks) F 19 
where 
Fy y? sin? (ry/h 
F, = y*? sin (ry fh) sin (2ryh 


satisfies the boundary conditions. It is 
differentiate Eq. (19 


atives 


. [ h \ ] TH \ 
F ] - | COS hb 
L a y/ J 2 ay 


necessary to 
and substitute the proper deriv 
into Eq. (1 


a F ?0 





FIG. 2 
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BUCKLING COEFFICIENT FOR PLATE SUBJECTED TO 
UNIFORM THRUST (FIG. la) 
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Derivatives of wat x = 0 become 


O-w ’ BY 4 k? rr? 
= = —F(y) 3 |2 + 

Oe" Feud a y" | 

O'w ‘ h\4 k? 3? k? xr? 
= Fy) ) 3 + 

ll z a y4 16 


(O7w Oy"), qo = 0? F(y) Oy? 


(0‘w/Oy*),-9 = O'F(y)/oy4 


( O*w ) (*) f O'F(y) 1 2 in kh? 1 OF(y) 4 2 k?3 7 if k?n? |) 
‘ . = <— : 7 2wr — © "(y) : 2+ 
Ox*Oy"/ 0 a { Oy ¥* } Oy y? ! 4" L 1 |! 
The derivatives of F(y) are written as 
7 wT ., say 4 BY 
F, _ y* sin + 2 y sin? 
h h 1 
= T , 2ry 4 2ry 5 €Y 
F, = 4 y sin + — ¥* cos + 2 sin? 
h h? f h 
eer 7” 2ry P . 2ry T 2ry 
F, = 94 cos — 32 y sin — vy? cos 
1 h 1 h h' h 
Re Samy 4 y? | ® in 4 Sr | Ory 
"9 = y]| cos — cos _~ sin sin 
} h 2 h h 1 h 
7 2ry . wy . omy Ty OmTy ys" Ty OOTY 
I ~ — sin + 3 sin + | cos — cos — — | cos — Jou — 
h r h h h 2h? h h 


67° 


pi anf ‘4 w [ . 
My = y sin 
¥ h | - fh Is 


If Eqs. (22) are substituted into Eqs. (21) and the resulti 
Pa + Oar aia 
Ra, a Sav — 


Ty 
— cos j + 9 cos 


h? a l 


where 
.- [k\*.... . 
P = %08:435 k?T — 42.052 + 39. 
a 
Q = 1736.350 — 12.5668 
i —_ 
R = 49.965 kT + 627.583 + 19.7 
a 
=i. 
S = —49.965 k?T — 1130.014 — 9% 
a 
and where 
8B = Nh?/D 
W = 2+ (wk?/4) 
T = 1+ 0.205622. 


Values of 8 for various ratios of h/a may be obtained 


from the condition 


ir 
R 3| 
Fig. 2 gives values of 8 for practical ranges of the param- 


Ty 
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ys" sm) 


» 7 

Ty Sry | 

* cos — S1 cos = | 
2h h H | 


: . ay] 4 
h h 


ng derivatives are substituted into Eq. (1), there results 
0, at y = h/2 


0, at y = 2h/3 


2.9348 


eter ha. In the case of an equilateral triangle 


B = 88.2 25 


(3b) CLAMPED PLATES UNDER SHEAR AND THRUS! 


For the load condition shown in Fig. 1(b) 

N, = —N(a/h)?, N, = WN 
If the first terms in brackets in Eqs. (24) are multiplied 
by —(a/h)*, the resulting equations will be sufficient or 





determining the critical value of NV. Fig. 3 gives values 


of 6 for different ratios of h ‘a. 
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FIG. 3 BUCKLING COEFFICIENT FOR PLATE SUBJECTED TO 
THRUST AND SHEAR (FIG. Ib) 


3) FINITE DIFFERENCE SOLUTION FOR CLAMPED PLATE 
UNDER THRUST 


In order to check the results presented in Part (3a) 
of this paper, the lowest critical buckling load of a 
uniformly compressed equilateral triangular plate is 
determined. Twenty-four interior points are selected. 
An equilateral triangular mesh is used.” 

The difference equations that result from writing 
Eq. (1) in difference form may be written matrically as 


ISOSCELES 
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(A + A\B)W = 0 265 
where 
A = matrix of coefficients of A?A°w terms 
B = matrix of coefficients of terms on right side of 


Eq. (1 
A = Nh? 50D 


Eq. (26) is in standard eigen-value form and the solu 
tion for \ may be obtained easily by using an Electronic 
Digital Computor. The value of ’ was found to be 
1.949. Therefore 

Ner = 97.45(D/h? (27 


This value is about 10 per cent higher than the value 
obtained in Eq. (25). 


(4) CONCLUDING REMARKS 


The curves presented in Figs. 2 and 3 should enable 
the engineer to obtain a rapid approximate answer for 
the lowest critical buckling load of isosceles triangular 
plates. It appears that the values of Ne, given by the 
curves are on the safe side in the practical range of the 
values of ha. 
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A Fraction Series Solution for Characteristic 
Values Useful in Some Problems of Airplane 
Dynamics 


E. L. LEPPERT, JR.* 
Lockheed . fircraft Corporation 


ABSTRACT 


Means for solving the determinantal equation D(A) =| A - 
kd | for its characteristic values is described. The method has 
been developed with particular regard for its use in solution of 
flutter problems and to its use with punch-card or high-speed 
electronic computing equipment. 

The roots are obtained by evaluation of a function consisting 
of a series of fractions and having the same roots as the poly- 
nomial of the The 
such a manner that its coefficients are readily calculated such 


determinant fraction series is formed in 


that the roots may be obtained to a high degree of accuracy 
A simple example illustrates the procedures 


D(X) = A * kn = Qy\ Rin ay 
(> (DD 
a> a 
Ani a 


having real or complex coefficients a@,, is required in the 
Wayland! de- 
scribes a number of methods for expansion into poly- 
The 


relative merits of the methods described are compared 


solution of many physical problems. 
nomial form preliminary to solving for the roots. 


by him on the basis of the number of operations re- 
quired. The important question of the requisite ac- 
curacy is not resolved. 

In the following development, a function consisting 
of a series of fractions? 


n+l 


n+l 


is utilized to determine all the roots to the same accu- 
racy. The X, are arbitrary values for substitution into, 
and subsequent evaluation of, D(A); the values of the 
C; coefficients are determined such that the roots of 
D(X) and f(A) are the same. 

The fraction series is the quotient of two polynomials. 
The numerator is the polynomial of the determinant; 
the denominator is a polynomial having as roots the 
arbitrary A;. Evaluation of the coefficients C; is made 

Received November 8, 1950. Revised and received July 8, 1954 
* Research Specialist 


9 
° 


y The method has been adapted to the IBM 701 Evaluation 
the 15 complex roots and their corresponding cigen-vectors for 
16 degree of freedom (one body mode without stiff comp] 
flutter determinant is accomplished in 59 min., including 
complete check. Evaluation of the roots without check requir 


29 min. 


INTRODUCTION 


S LUTION FOR THE eigen-values of the determinantal 


equation 


143 ay 

sic Ror a a 
a — ksr a = 
a a, = kX 


according to Peirce.* In general form, 


D(X;) 
n+ 1 9 

IT (A; — A 
z ] 

For the special case where the substitution ), is a 
root, so that II(\, A;) 0: PIX.) = 0 and €. =8 
The function f(A) then has the same roots as the deter- 
minant out. While selection 
of the \; may be completely arbitrary, evaluation of the 


with one root divided 
series to determine the roots is improved if the substi- 
tution values approximate the roots as closely as 
possible. In the absence of any knowledge of the 
distribution of roots, selection of a geometric progression 
has been found advantageous for some vibration prob- 


lems. 
EXAMPLE 


An example illustrates the simplicity of the opet- 
A determinant having the eigen-values 0.2, 
The polynomial to which such 


ations. 
1.7, and 7.0 is assumed. 
determinant expands is 


Q 


D(A) = (A — 0.2) (A — 1.7) (A — 7.0) 0 


Presuming no information other than the order ol 
the polynomial, the arbitrary values \ are chosen 
0,1, 2,4, and 8. Evaluation of the coefficients for the 
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SOME PROBLEMS OF 


_—" 
TABLE | 
—T - : —— i esdcmcemes = c/ 
a \ j -Ai) | D(A j) J 
2 4 ~ ij _DIAi) _ | 
J | 0 Aj -Ai) 
+ + + + + = 
-3 -7 | -2 3.36 -.1600 | 
3 2 5 i2 2.70 2250 | 
| 
4 2 0 4 24 26.22 0925 
| 
g 7 6 + 0 168 49.14 2925 | 
: -238 | 
: - l | — 





fraction series 1s shown in Table 1 from which 

0.2250 

A -—2 
1.0925 


A-—-4 A-—s 


0.1600 


A— | 


Examination of the values of the determinant D()\) 
reveals changes in sign between O and 1, between | and 
2 and between + and 8. Three methods for obtaining 
the roots are shown 

1) Substitution of assumed values of \,, calculation 
of the corresponding values of /(A) (see Table 2) and 
plotting the resultant curve as in Fig. | until a crossing 
is obtained is rapid and well suited to desk calculation. 
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(2) An iterative solution, wherein an assumed )\,4 is 
substituted in all the fractions but one and a value i, 
obtained for that one such that /(A 0, is often pos 
sible. Substitution of Av into the fractions in place 
of A4 and solution for a new Ace generally results in 
rapid convergence to a root lable 3 lists values of 
Ac for a number of A, to illustrate rapidity of con 
vergence toward the root from widely separated as 
sumed values. 

(3) A modified iteration procedure requiring calcu 
lation of the Ac for two trial values of A, in the vicinity 
of the root and assuming linearity of the function has 
Postel for use by the IBM 


Card Programmed Calculator to solve for the eigen 


been developed by E. E. 


values of the complex flutter determinant. Corre 
sponding to the two trial values, two calculated , 
values, C; and Cs, are obtained. A third trial value, 


3, is then obtained with the linearity formula 
Co — dre 
Cc, -— C, 
Ar — A 


] - 
and a corresponding C3; is calculated. The linearity 
formula is then again applied to obtain a fourth trial 
for which the corresponding C; is calculated. The 
procedure is repeated until sufficient convergence of 
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TABLE 2 

AA (1) 2 (3) 4 f(A 
0 +0.1600 +0.1125 0.2731 0.0366 0.0372 
0.1 +0.1778 +O.1184 0.2801 0.0369 0.0208 
0.5 +0.3200 +0.1500 0.3150 0.0390 +0.1160 
l € t 
l+e 
1.3 0.5333 +0.3214 0.4046 0.0437 0.6516 
1.5 0.3200 + 0.4500 0.4370 0.0450 0.3520 
1.8 0.2000 + 1.1250 0.4960 0.0472 +0.3818 
2 € 
2 € 
2.2 0.1333 1.1250 0.6100 0.0505 1.9188 
2.4 0.1140 0.5625 0.6850 0.0523 1.4138 
3.0 0.0800 0.2250 1.0925 0.585 1.4560 
3.7 0.0593 0.1325 3.6600 0.0681 3.9199 
3.8 0.0572 0.1250 5.4800 0.0696 5.7318 
i—e 
{+e 
6.0 0.0320 0.0750 +0.5480 0.1461 +0. 2949 
6.5 0.0291 0.0500 +0.4380 0.1950 +0.1639 
7.5 0.0246 0.0409 +0.3130 0.5850 0.3375 
7.25 0.0256 0.0429 + 0.3370 0.3900 0.1215 


r, and C,, is obtained. An example utilizing numeri- 


cal values from Table 3 follows: 


Ai = 1.30, C, = 1.77 
Ae = 1.50, C2. = 1.718 
As = 1.673, Cs; = 1.7014 4, = 1.699 


DISCUSSION 

The procedure is used at Lockheed as one of several 
methods of solution of the flutter determinant. In 
many cases, all the eigen-values are obtained by it 
with the same amount of computation as is required 
to obtain two or three roots by other methods. Utiliza- 
tion of the latest available IBM equipment to evalu- 
ate the determinants is a basic premise of the method. 

An important feature of the method is the accuracy 
of the solutions. Each coefficient of the fraction series 
is obtained to the same accuracy as the evaluation of 
the determinant. They are not obtainde by differ- 
ences of large numbers, a condition that requires com- 
putation with an exorbitant number of significant fig- 
ures in other methods. The resultant eigen-values, 
whether large or small, are thus obtainable to the same 
accuracy as the value of the determinant. 

Solution for the eigen-values of all real determinants 
of 16th order, and of complex flutter determinants of 
Sth order, is readily accomplished with presently avail- 
able punched card equipment. Routine evaluation 
of each Sth order complex determinant requires | hour 
for each starring reduction on the IBM 604. From 
1.8 to 2.0 hours on the Card Programmed Calculator 
are needed to form the fraction series and to iterate to 
all the roots. 

The Lockheed Mathematical Analysis Department 
has programmed the method for the IBM 701 ‘“De- 
fense Calculator.”’ Evaluation of an Sth order com- 
plex determinant requires 15 sec. In an application to 
16th order complex flutter determinants, the 15 eigen- 
values and their corresponding eigen-vectors are ob- 
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TABLE 3 
AA r 
0.1600 
0 + 9.1125 0.273125 0.0365625 0 LO 1 RR, 
IY l = 
0.1600 
0.5 + 0.1500 0.31214 0.0390 0 0% 
AC 1 Ps 
0.1600 
0.1 + 0.11842 0.280128 0.037025 
ny, ] 
0.1600 
0.2015 + 0.1251 0.2876 0.0375 0 0.2 
XC l 
0.2250 
1.3 0.5333 +4 0.4046 0.0437 0 177 
x 2 
0.2250 
1.5 0.3200 +4 0.4370 0.0450 0 718 
nu 2 
0.2250 
1.8 0.2000 + 0.4966 0.0472 0 607 
x 2 
0.2250 
1.673 0.2377 + 0.4695 CG .0462 0 ™ 
ny, 2 
7 0.2286 + 0.7500 0.4750 0.0464 0 


tained in an average elapsed time of 59 min. for each set 
The modes include one body mode with no stiffness 
The time includes a complete check. Evaluation of 
the roots without check requires 29 min. 

Improvements on the 701, which will be available on 
production within a year, will allow evaluation of the 
larger degree of freedom problems in about one quarter 
the time. 

The method has been utilized to explore the effect of 
variation of tip tank parameters on all of the wing 
flutter modes of a large military airplane. By utilizing 
primitive modes having motion of the tip tank in only 
one twist and one bending degree of freedom, it is 
possible to economically determine values of the deter- 
minants for a wide range of tip-tank conditions, utiliz- 
ing a single set of assumed \. Formation of the corre- 
sponding fraction series and solution for its eigen- 
values thus represents almost all the additional work 
required to obtain a complete set of solutions for each 


tip-tank condition. 


CONCLUSIONS 


A means for solving the determinantal equation for 
all its eigen-values is described. The method is of 
particular value in solution of flutter problems with 
punch-card equipment. A particular feature is its 
ability to obtain all roots with good accuracy. It has 
been found particularly useful in the study of the effect 
of variation of parameters in a problem having a large 
number of degrees of freedom. 
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Stability of Compressible Laminar Boundary 
Laver with Internal Heat Sources or Sinks 


GEORGE M. LOW* 
Lewis blight Propulsion Laboratory, NACA 


SUMMARY 


of the equations of the compressible ] uminar bound 


The 


4 solution 


laver with heat sources or heat sinks is presented 


¢ temperature and constant Prandtl Number. It is assumed 

t the additional mass associated with the heat sources or sinks 

is negligible 

Stabilitv calculations indicate that the boundary layer can be 
stabilized by heat withdrawal near the aircraft surface (e.g., by 

yaporation of liquid droplets) or by heat addition near the outer 
ige of the boundary layer (e.g., by condensation of a vapor 
In addition, heat withdrawal near the surface is shown to pro 


vide an effective means of surface cooling 


INTRODUCTION 


| aes AERODYNAMIC HEATING of high-speed aircraft 
and missiles leads to many new design problems. 
If flight speeds continue to increase at the present rate, 
the steady-state skin temperatures encountered will 

the metals 
used in aircraft construction (see Fig. 1). 
of these metals is to be continued, adequate means of 


for 
If the use 


soon exceed allowable temperatures 


cooling will have to be provided. 

Coexistent with the cooling problem is the problem 
of boundary-layer transition. If the boundary layer 
becomes turbulent, approximately ten times the amount 
of coolant is required to maintain a given surface tem- 
perature than that needed if the boundary layer re- 
As yet, little is known about the mech- 
flow; 


mains laminar. 


transition from laminar to turbulent 
presumably, if the turbulence level is low (as it will be 


in flight), transition is a result of the instability of the 


anism of 


laminar boundary layer to small disturbances. 

The stability of the laminar boundary layer in a 
compressible fluid has received a great deal of attention 
inrecent years. In particular, Lees and Lin! and Lees? 
have established criteria for the determination of sta- 
bility limits. In these papers it is suggested that, if 
sufficient heat is withdrawn at the surface, the bound- 
ary layer can be stabilized for all Reynolds Numbers. 
Detailed calculations by Van Driest* show that com- 
plete stability to two-dimensional disturbances can be 
Mach The 
amount of cooling required for stability was found to 


attained at Numbers between 1 and 9. 


be of the same order of magnitude as required to main- 


tain the wall temperature at a reasonable level. These 


> 4 “ . ° . . . 
Presented at the Aerodynamics Session, National Summer 


Meeting, IAS, June 21-24, 1954, Los Angeles 
* Head, Mechanics 
Division 


Fluid Section, Supersonic Propulsicn 
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results have been at least qualitatively substantiated 
by Sternberg,’ who observed laminar boundary layers 
at Reynolds Numbers as high as 50 X 10°. 

In the past, two methods of cooling, and thereby of 
stabilizing the laminar boundary laver, have been sug- 
gested: surface cooling, accomplished by passing a 
coolant along the inside of the skin, and transpiration 
cooling, accomplished by passing the coolant through 
the skin. In the present report a third method will be 
investigated —that of cooling by heat sinks within the 
boundary layer. the 
boundary layer by use of heat sources will also be in- 


The possibility of stabilizing 


vestigated. Physically, the evaporation of liquid drop 
lets in the boundary layer will act as a heat sink; 
whereas condensation of vapor or burning of fuel will 
be a source of heat. 

A prerequisite for the calculation of stability limits 
is a solution of the applicable boundary-layer equa 
tions. These equations are therefore derived and solved 
for several particular distributions of heat sinks or 
sources. The Lees-Lin stability criteria are then ap 


plied to these solutions. 


SYMBOLS 


The following symbols are used in this report 





Fic. 1 





ayx** = polynomial describing wall temperature distribu 
tion (ayx** a ayx* + aox*? 

Cc = constant of proportionality in viscosity temper 
ature relation 

C; = average friction drag coefficient 

& = the factor ¥ t,/t,, (t re le t+ S), a function 
of x 

U = phase velocity of disturbance, divided by u 

Cc = local skin-friction coefficient 

C; = specific heat at constant pressure 

f = solution of momentum equation 

g = function describing temperature deviation due to 
heat sources (Eq. 18 

u 3000 
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= thermal conductivity 


9 = reference length 

| = Mach Number 

P) = Prandtl Number, uc, /k 

) = static pressure 

QO = heat addition per unit mass per unit time by source 


(or withdrawal by sink) 
(?" = OL Cpl ul 
= local rate of heat transfer 


R = gas constant 

Re = Reynolds Number, p..xu, /m 

, = temperature function (see reference 5 
S = Sutherland’s constant (S = 216°R 


t = static temperature 


t/t 
ul = velocity in x-direction 
it = Uu/u 


= velocity in y-direction 
= coolant flow rate 
distance along surface measured from leading edge 


\ = x/L 


iy = temperature function (see reference 5 
y = distance from surface measured perpendicular to 
surface 

+ = ratio of specific heats 

n = characteristic variable defined by Eq. (9 
= variable defining heat source location 

im = coefficient of viscosity 

rm = p/p 

v = kinematic viscosity, u/p 

p = mass density 

T = shearing stress 


é = local conditions at outer edge of boundary layer 
lw = insulated wall 
= conditions at wall or surface 
= conditions at x = 0, y = 
x, ¥,7 = partial differentiation with respect to x, y, or 
Superscripts 
* = dimensionless quantities 


= total differentiation with respect to 
A bar over a quantity indicates integration from zero to 7 


BOUNDARY-LAYER ANALYSIS 


Differential Equations and Boundary Conditions. 
The steady laminar flow of a viscous compressible 
fluid in a boundary layer on a flat plate is governed by 
the momentum equation 


uu, + vu, = (1/p) (ut,), Cn) 


the equation of continuity 


to 


(pu), + (pv), = O ( 
the energy equation 


Cpt, + Cyvt, = (1/p) (Rt), + v(u,)? + Q (3) 
and the equation of state 
pRt = p = constant (4) 


The last term in the energy equation describes the heat 
addition (or withdrawal) per unit mass per unit time. 
These equations are subject to the following boundary 


conditions: 
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u(x,V0) =0 u(x, ©) =u 
v(x, 0) = 0 he, @) = fF 
t(x,O) = 1t,(x%) (for a noninsulated wall 
or 
Eis, 0) = 0 (for an insulated wall 
Assumption and Limitations.—The following assum 


tions and limitations are imposed in order to simplif 
the solution of Eqs. (1) to (4 

(1) The viscosity and temperature are linearly y 
lated as follows: ; 


a = ee ( 
(Asterisks denote dimensionless quantities as defined 
under the section, Symbols.) Chapman and Rubesin 
have shown that solutions of the boundary-layer 
equations based on this viscosity-temperature lay 
agree well with more exact solutions at Mach Numbers 
up to about 5 if the constant C is determined by match 
ing Eq. (6) with Sutherland's relation at the solid 


boundary: 


This assumption should also be reasonable for flows with 
heat sources, provided the resulting temperature devi 
ation is not large. 

(2) The additional mass associated with the heat 
sources has a negligible effect on the boundary-layer 
development. The validity of this assumption is dem- 
onstrated a posteriori. 

(3) The surface temperature can be represented by a 


polynomial in x: 
‘.* = f ia + anx*® \ 


where /;,,.* is the dimensionless wall temperature for 
zero heat transfer and the repeated index .V indicates 
a summation. 

Solution of Equations._-The present problem and the 
problem solved by Chapman and Rubesin® differ onl 
in the last term in the energy equation. It is therefor 
only natural that the solution of the equations pre 
sented herein follows along the same line as given 1! 
reference 5. Priefly, the procedure is as follows 
The system of four equations [Eqs. (1) to (4) ] is reduced 
to two equations by absorption of the continuity equa- 
tion and the equation of state into the other two equa- 
tions. The dependent variables in the remaining 
equations are a newly defined stream function and 
temperature. Further reduction of these equations 
is obtained by introducing new independent variables 


x* and n, defined as follows: 


1 | yy , 
= | ly 
’ 2 y p G J, i” - 


The dependent variables in this coordinate system are 


fand ¢ and are defined by 
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(10) 


[he resulting equations and associated boundary con- 


ditions aré 


+4 Prit,* — 2Prf’x*le* = 
7 . 2 " 
= M_2Pr(f")2 — 4Prx*O* 
| 
2) 
fixe". oO = f + dyx*N 
“(x*, =1+ | nts”, x dx* 
J 0 


The momentum equation [Eq. (11)] is, of course, the 
well-known Blasius equation. The energy equation 
Eq. (12)] is, with the exception of the last term on the 
right-hand side, identical to the energy equation solved 
by Chapman and Rubesin. Pecause the equation 1s 
linear, only a particular solution satisfying Eq. (12) 
need be added to the temperature function of reference 
5 in order to obtain the complete solution. 

Up to this point the derivation has been quite gen- 
eral. Before a solution of the energy equation can be 
obtained, however, a particular form for the heat-source 
At this time, there- 


fore, the assumption is made that Q* is a function of 


distribution need be assumed. 


n alone. In order to solve Eq. (12), it is divided into 


two linearly related parts, 

t* = t,* + t.* 3) 
where /,;* and /:* satisfy the following equations: 
ho» + Prft,,* — 


2Prf'x*te* = 


ek ” 
- , M .?Pr(f 


(14) 
t*(x*, O) = tio*® + ayx*® — to*(x, 0) 
i. *(x*, )= ] 
fon.* + Prfts,* 2Prf'x*te.0* = —4Prx*O* 
t ¥(x* UU) = Gg (15 


te* (2*. eo) = x*0*| 20 ) 


Eq. (14) therefore describes the temperature distribu 
tion caused by wall heat transfer and by frictional heat- 
ing. Eq. (15) describes the temperature increment due 
to heat sources or sinks. The wall heat transfer in 
Eq. (15) can be taken as zero without loss of generality, 
because arbitrary heat-transfer conditions are included 
in Eq. (14). The solution of Eq. (14), as given in refer- 


ence 5, 1s 


ere l ' 
Le*,s) = 1+- M ?r(n) + ayx**¥Vy(n) (16) 
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The solution to Eq. (15) is found to be 
to* = 


1PrQy*x*g(n 17 


where g is the solution of the ordinary differential equa- 


tion 
g” + Prfg’ 2Prf'g VY" Vo" 
IS 
o'(0) = 0, gl = 0 1PrQ,* 


and where Q* is a measure of the total source or sink 


strength, defined by 

Q,* = { OG" d y 

‘ Jo © an If 
The functions f, r, and }y and, where required, their 
derivatives and integrals, are presented graphically in 
reference 5. When more accuracy is desired, tabulated 
values of f and r can be found in reference 6, while 
numerical values of }’y can be obtained from the fol- 


lowing relations: 
Hyiln . 20 


where the functions on the right-hand side are tabu- 


lated in reference 6. The function g(n) was obtained 
on an electromechanical differential analyzer for nine 
specified distributions Q*(n). These distributions can 
be classified as ‘‘broad-band’’ distributions, defined by 


OF =10<n< 7 


Oly )» 


Sketches of these distributions appear as inserts in 
ey 
Fig. 2. Values of g, g’, and g { g dn are presented 
) 
(Pr = 0.72. 


Boundary-Layer Characteristics. 


in Figs. 2, 3, and 4. 
Fecause the bound 
ary-layer characteristics are obtained in the same 
manner as shown in reference 5, only final results will 
be presented herein. 

Velocity and temperature profiles are obtained as 


functions of n by 


u* = (1/2)f'(m) 23) 
y-— 1 
*=%4*+h* = 1+ - . M .? r(n) + 
ayx** Vy(n) + 4PrOQo*x*e(n) (24) 


where the coefficients ay are obtained by fitting a poly- 
nomial to the wall temperature distribution 


Anx** = ty* — tin* 25 


The insulated wall temperature is 


M2 r(0) + 4PrQo*x*g(0 26) 


a 
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(b)NARROW-BAND DISTRIBUTIONS OF SOURCES OR SINKS 


Fic. 2. Variation of temperature increment due to heat sources 
or sinks 


If it is desired to obtain profiles in the physical x-, y- 
coordinate system, then the following relation between 
yvand n is required : 

M ??(n) + 


ayv*®* Vy(n) + iProuts*acn) | (27) 


(A bar over a quantity denotes integration from zero 


ton.) The wall shearing stress is 
On f" (0) ; lo Uw 
To (. ) =~ a oe : (28) 
OV) x | tx 
where C,(x*) = uy*/t,* is equal to C for the special 


case of constant wall temperature. Local and average 


skin-friction coefficients are defined as follows: 


T moO) C¢C, 


Cr, = = (29) 


(1/2) 6.0.” 2 VCRe 


| 
Ca = ‘a | Tx dx 
(i/o) 0.x 70 
which, for constant wall temperature, becomes 
Cr = 1.328+>/C/VRe (30) 


Thus it can be seen that heat sources affect skin 
friction only indirectly. In general, heating will alter 
the wall temperature, which in turn affects the constant 
C. The resulting change in skin friction is expected to 
be small. 
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The local rate of heat transfer is found from the tem. 
perature gradient by the following expression 


q —kiy(ty)w 


The effect of heat sources or sinks on heat transfer js 
expected to be appreciable, because they have a first 
* 


order effect on ¢;,,*, which is directly related to ayx* 
For example, it is expected that the transient heating oj 
a missile could be reduced appreciably by the use of 


heat sinks in the boundary layer. 


UsE OF HEAT SINKS TO COOL THE SURFACE 


The use of a broad-band distribution of heat sinks to 
cool the surface was investigated. In particular, the 
amount of heat withdrawn by heat sinks over an in 
sulated surface [Fig. 5(b)] was compared with the 
amount of heat withdrawn at a wall [Fig. 5(a)] with 
the same surface temperature. The relative heat- 
transfer rates are shown in Fig. 6. The significance of 
these rates of heat transfer can be interpreted in terms 
of two alternate methods of cooling. If, for example, 
water is used as a coolant and if it is fully evaporated 
in both systems, then the relative coolant flow rates 
will be the same as the heat-transfer rates; under these 
conditions, wall cooling is slightly more efficient than 
cooling with heat sinks. 

On the other hand, conditions may be such that 
evaporation can take place within the boundary layer, 
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whereas Within the wall the coolant cannot boil. 
When compared on the basis of coolant weight re- 
quired, the heat sinks then provide a much more effi- 
cient cooling system than that provided by wall cooling 


Fig. 7 


STABILITY ANALYSIS 


The stability of the laminar boundary layer on a flat 
plate in a compressible fluid has been investigated by 
Lees and Lin. As for the incompressible case, the 
stability differential equations are obtained by linear- 
izing the equations of motion and energy for a flow with 
time-dependent fluctuating components superposed on 
the steady mean-flow quantities. These equations, in 
general, lead to a characteristic value problem relating 
the wave velocity and wave length of a disturbance to 
the parameters that govern the shape of the velocity 
and temperature profiles. 

There is, however, one difference between the incom- 
pressible the stability 
In the compressible case, distinction must be made 
By 


and compressible problem. 
between supersonic and subsonic disturbances. 
definition, the x-component of the phase velocity of a 
supersonic disturbance relative to the free-stream ve- 
locity is greater than the free-stream speed of sound, 
and conversely. When a disturbance is supersonic, the 
characteristic value problem does not arise, and Lees 
and Lin have reasoned that amplified supersonic dis- 
turbances cannot exist. It therefore follows that only 
subsonic disturbances—that is, disturbances for which 
the dimensionless phase velocity c is greater than 1 — 

M,)—are relevant in boundary-layer stability cal- 
culations. 

Lees* suggested the possibility that, if sufficient heat 
is withdrawn from a supersonic laminar boundary layer, 


the flow will be stable at all Reynolds Numbers. This 
(y — 1) M27’ + 

— } a 
1+ [(y — 1)/2)M 


The following procedure was employed in the deter- 
mination of stability limits: For a given wall temper- 
ature distribution, a value of Qy* was found such that 


the left-hand side of inequality (35) vanishes when 
c= 1 — (1/4,). 
determine whether the criterion was satisfied when 
c>1— (1/M)). 
was the case. In the few cases when the criterion was 
not satisfied, Q)* was adjusted by trial and error until 


» 
oe 
This value of Qo* was then used to 


In general, it was found that this 


the correct value was determined. 

The basic stability criteria, Eqs. (32) and (33), were 
derived for the flow without heat sources or sinks. It 
has been shown in the past’ that stability criteria de- 
rived for flat-plate flows apply equally well to flows with 
pressure gradients, provided the local velocity and tem- 
perature profiles are considered. In other words, pre 
vious history of the flow is important only insofar as it 


allects the profile shapes. By analogy, it is assumed 
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means that no self-excited disturbances can exist 


which satisfy the stability differential equations, their 
boundary conditions, and the requirement that c be 
M, 


“complete stability” 


greater than 1 — (1 In the present report only 


conditions for this are investi 


gated. 

Solutions of the stability problem fall in two general 
categories—inviscid solutions and viscous solutions. 
In the inviscid solutions, as the name implies, the effect 
of viscosity is neglected; these solutions apply only at 
In the 


the effect of viscosity is considered to the first order; 


infinite Reynolds Numbers. viscous solutions, 
these solutions therefore apply at large, but finite, 
Reynolds Numbers. In general, the inviscid solutions 
are much easier to obtain than the viscous solutions. 
However, care must be taken in interpreting inviscid 
results. A flow that is stable at Re will usually 
not be stable at all Reynolds Numbers. It was found 
in reference 3 that, as far as complete stability is con 
cerned, the viscous and inviscid solutions vield similar 


In 


the present report, therefore, only the inviscid solution 


results at Mach Numbers greater than about 1.5. 


is considered. T 
The inviscid stability is determined by the shape of 


the velocity and temperature profiles in the boundary 


layer. In particular, amplified disturbances will not 
exist if 

(0/Ov) [p(Ou Oy)] < 0 32 
when 

c=u*>1 — (1/M, (33 


In terms of the coordinates used in the present report, 
Eq. (52) becomes 


(mn? a"). = (EF/ FP) <0 34 


or, by use of Eqs. (23), (24), and (11 


, 


*\ V5 , + SPrQo*x*g 


anx** YY. + 1PrQy*x*¢g 


herein that the stability equations derived for flows 
without heat sources apply also to the flows treated in 
the present analysis if the proper velocity and tem 


perature profiles are used. 


It should also be noted that instability to two 


dimensional disturbances only is considered. It was 
shown by Squire® that, for incompressible flow, two- 
dimensional disturbances are always more destabilizing 
than three-dimensional disturbances. Dunn and Lin’ 
have recently shown that Squire’s theorem cannot be 
extended to sufficiently high 
speeds, the boundary layer can never be completely 
with to all 

However it has been shown experimen 


supersonic flows; at 


stabilized respect three-dimensional dis 


turbances. 
use of the viscous 


10(a) was calculated by 


This point is seen to yield essentially the same result 


t One point in Fig 
solution 
as the corresponding inviscid result 
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it sources; g 


tally? that surface cooling has a large effect on sta 
bility. It is believed, therefore, that the present re- 
sults at least qualitatively describe the effect of heat 
sources and sinks on boundary-layer stability. 


EFFECT OF HEAT SOURCES AND SINKS ON STABILITY 


rhe stability of the laminar boundary layer is gov- 
erned to a large extent by the distribution of angular 
momentum throughout the layer. The theory states 
that, if the angular momentum increases with y in that 
portion of the boundary layer where u* > 1 — (1 /J/,), 
disturbances may be amplified. At supersonic speeds 
such an increase in p(Ou Oy) can be caused by the 
large negative temperature gradient in the central part 
| the boundary layer [see Eq. (34)]. If the tem- 
perature gradient can be kept small enough so that 
at > =f (36) 
the flow will be stable. This distortion of the tem- 
perature profile could be accomplished by withdrawing 
heat at or near the wall or by adding heat at a distance 
irom the wall. 
In Fig. 8 are shown three temperature profiles that 
will just stabilize the boundary layer at a Mach Num- 
ber of 3.69. Also shown in this figure is the profile for 


the unstable flow over an insulated surface. It is ap- 
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parent that whether the stabilizing is done by wall 
cooling, by heat sinks, or by heat sources, the effect 
is always the same: The central part of the stable 
profiles is less steep than the central part of the un 
stable profiles. The outer portion of the profile with 
heat addition is fairly steep, but there the value of f is 
also large, so that condition (36) 1s not violated 
It was stated in the section entitled ‘‘Assumptions 
and Limitations’ that the additional mass associated 
with heat sources or sinks is negligible. The mass 
flow of water or steam required to stabilize completely 
the boundary layer is shown in Fig. 9 as a percentage 
of the total boundary-layer mass flow. The addition of 
mass flow ranges between 0.01 and 2.00 per cent. By 
assuming that the additional mass is diffused like the 
temperature, it was found that the maximum percentage 
change in u” u’ is about three times as large as the per 
centage change in mass flow. Thus the effect of neg 
lecting mass addition appears to be small. 

Stability with Heat Sources and Sinks.—In Fig. 10 
with heat sources or sinks on an insulated plate as com 
pared with stability by wall cooling are presented. 
The broad-band distributions of heat sinks |Fig. 10(a) } 
are effective only if they are at or near the wall; at a 
distance from the wall the heat-transfer rate becomes 
prohibitive. Narrow-band distributions of heat sinks 
are effective only when they are extremely near the 
wall [Fig. 10(b)]. When the narrow-band distribu 
tions are placed at a distance from the surface, then heat 
addition is required, as indicated by the negative sign of 
the quantity ¢ qo. 

When the narrow-band distributions are placed in a 
critical layer in the central portion of the boundary 
layer, an infinite amount of heat addition or with 
drawal is required to stabilize the boundary layer 
The location of this critical layer will, in general, de 
pend on the Mach Number and the exact shape of the 
distribution. From Eq. (35) it is seen that an infinite 
amount of heating or cooling is required for an insulated 


surface when 


2g’ + gf = 0 (37) 
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at that point in the boundary layer where 
u”™ = 1 — (1/M,) 


If Eq. (37) is solved for g and the result is substituted 
into Eq. (18), a distribution is obtained which cannot 
stabilize the boundary layer at any Mach Number. 
This distribution is shown in Fig. 11. 

The relative efficiency of the stabilizing systems de- 
When 


water is used as a coolant, heat sinks can be much more 


pends upon the method of cooling or heating. 


effective than wall cooling if evaporation can only be 
realized within the boundary layer but not beneath 
the skin. If, on the other hand, evaporation can take 
place both within the wall and in the boundary layer, 
then wall cooling appears to be the more efficient 
method. Another possible system of stabilizing in- 
volves the simultaneous use of wall cooling and heat 


addition in the boundary layer [Fig. 5(c)]. The fol- 
lowing example of this system was calculated: Water 


at 40°F. and at sea-level pressure is passed along the in- 
side of the skin, where it is evaporated and the resulting 
steam is raised to the surface temperature. The steam 
is then ejected in the outer portion of the boundary 
layer (ym) = 2.5), where it reaches the local velocity by 
a constant total-enthalpy process. Finally, condensa- 
tion of the steam takes place. A ratio of w/w» was 
calculated, where w is the coolant flow rate for this sys- 
tem and where wy is the amount of water required to 
stabilize the boundary layer through evaporation within 
the wall. The following results were found for an alti- 


tude of 50,006 ft. : 


M WwW ral 
.. 3 0.90 
es: 0.74 


At high Mach Numbers, therefore, a sizable reduc- 
tion in coolant weight can be realized by reusing the 
evaporated coolant as a source of heat in the outer por- 
tion of the boundary layer. 

When the boundary layer is stabilized by wall cool- 
ing, complete stability cannot be attained at Mach 
Numbers greater than about nine.* This limitation 
does not exist when the boundary layer is stabilized by 
use of heat sinks, although the amount of heat that 
must be withdrawn at high Mach Numbers becomes ex- 
tremely large. For example, at 1/7, = 13, eight times 
the amount of heat removal is required for stability than 
that required when 1/7, = 6.5. This result, of course, 
is only approximate, because neither the stability 
theory nor the boundary-layer equations used herein 
apply at these high Mach Numbers. 


PRACTICAL ASPECTS 


In all calculations it was assumed that the heat 
sources or sinks are evenly distributed in the direction 
of the flow and that their strength varies essentially as 
Vx. In practice it would probably be desirable to 
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inject the sources or sinks at one, or several, distinc 
stations and to obtain the proper distribution throug} 
the time lag inherent in the process of evaporation o; 
condensation. It is not known how closely the the, 
retical problem duplicates the practical conditions, 
he question naturally arises as to whether it 
feasible to inject the heat sources or sinks without thy 
injecting mechanism itself being the cause for trangj 
tion. It appears possible that, if the injecting mech 
anism is placed in a stable part of the boundary layer 
and if it does not generate a large wake, all disturb 
ances originating there will be damped. If this possi 
bility is confirmed experimentally, the relative ad- 
vantage of heat sinks for cooling purposes and of sinks 
or sources for stability purposes will depend primarily 
on the possible weight-savings. . 
A comparison on a weight basis with other systems 
was not undertaken, because it is difficult to define a 
system without having a specific airplane in mind 
Circulation of fuel past the surface, circulation of an 
expendable coolant past the surface, and transpiration 
cooling have all been suggested in the past. The 
aforementioned systems would be rather complex, be- 
cause the airplane skin would have to be either a com- 
plicated heat exchanger for purposes of surface cooling, 
or porous, for transpiration cooling. The use of heat | 
sinks might therefore represent a substantial weight- 
saving, because neither a double shell structure nor a 
porous structure would be required. 
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Note on the Unsteady Heat Regenerator 
Meyer and T. A. d’Ews Thomson 
Department of Aeronautical Engineering, University Sydney, 
Sydney, Australia 


December 7, 1954 


N MOST SUPERSONIC TUNNELS, control of stagnation pressure 
posed for the control of Reynolds Number, but in blow 
down tunnels of large working section and small running time, 
ontrol of stagnation temperature is similarly important The 
provision of such control by passing the air, upstream of the 
nozzle, through porous material has been discussed by Pillow, 


issuming that (1) the material is packed uniformly into a cylin 


irical duct with heat-insulated walls, (2 


heat conduction in the 
tir and in the material is negligible compared with the heat 
transfer between the two, and (3) the density is uniform through 
ut. Pillow shows that with finely constituted material a satis 

tory approximation for the purpose of stagnation temperature 
outrol is obtained by neglecting the difference between the tem 
eratures of the vir and the material at any fixed station, x, and 
time, f 

On the other hand, since the stagnation pressure needs to be 
similarly controlled, it is in general not the density, but the pres- 
sure, Which is constant in the regenerator, except for a small 
pressure drop through the porous material, which will be neg- 


ected in the following. The equation of continuity is 


Op/Ot + A pu)/Ox = 0 (1 


The heat bal 
ince for the segment of the duct between x and x + dx is 


where p is the air density and u the air velocity 


(hs + pc,)0T/O¢ + c,O0( puT)/dox = 0 (2 


where 7 is the temperature, /: is the density of the porous material, 
ind c, and s are the specific heats of air and material, respec 
With the pressure constant, 0(p7°)/0t = 0 and from Eqs 
1) and (2), 


tively 
O[(pu)?/(’ + 2/T)]/dox = O 


where X = hsR (pe, 
Stant 


, P being the pressure and R the gas con 
In general, the mass-flow rate m at the outlet will be 
constant and known, and, if A is the area of cross section and 
§\t) the value at the outlet of 


(pu)? = (m/A )0/O. 


ind from Eq. (2 


(p/R)0O/Ot + (m/A) (0/0 00/dox = 0 


or in terms of @ and ¢ as independent variables, 


p/R)dx/dt = (m/A) (062) (3) 


if .s i i i 
Here Ox /O¢ is the local rate at which any given temperature is 
transmitted through the duct. 

For stagnation temperature control it is desirable that the initial 


‘mperature distribution, 7)(x), of the porous material be uni 


nN the aeronautica Clence 


JOURNAL are presented in this special department. Publication 


expre 


The Editorial Committee does not hold itself 
ed by the correspondent 
form Phen the solution of Eq. (3) is 7(x, ¢ 7) as long as 
t < xt /l, where / is the length of the duct and 
tx = pO&Al/(Rm) = hAl [(s/e,) + (2p0/h)] /m 4 


rhis is the duration of the run at constant outlet temperature, 
and it is longer by a factor (s/c, + 2p0/h)/(s/Cp + po/h) than that 


for incompressible flow. Since 6(¢) = 6 during this period, 


Eqs. (3) and (4) give the equation 


} 


x/l = (0/6)'2[1 — t(T)/te] 


for the temperature distribution in the porous material at time 
tx, where ¢,;(7) is the time at which the temperature 7 occurs at 
the inlet. From this the regeneration required between runs can 
be estimated 

After this initial period, the rate of transmission of temperature 
decreases not only with the temperature but also falls during 


Integration of Eq. (3 
7 
? dr 


n(T) 97) 


the transmission at constant @ gives 


tx 0 = 
e 
by Eq. (4), where f(7°) is the time at which the temperature 7 
occurs at the outlet; therefore, 


dto/d@ = [0/00(t,)] (*dt;/d0 + (1/2)ts /O 
For tf; S tx, which includes the period of interest for wind-tunnel 
operation, 0(t;) = 9, and, hence, 
a. ; 
to — ty = te + 1/ate(0/0 — 1) 4 f (0;/8o — l1)dt 
0 


where 6,(¢) is the value of @ at the inlet. This shows how the 
delay in the transmission of temperature from inlet to outlet in 


creases as the temperature at the inlet falls 
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A Procedure for the Approximate Analysis of 
Buckled Plates 


Sara R. Boley 


Department of Civil Engineering and Engineering Mechanics, 
Columbia University, New York, N.Y 


December 8, 1954 


A METHOD FOR OBTAINING approximate solutions to the von 
Karman plate equations,' and in particular the effective 
width of buckled plates is presented below. It is an extension 
of a method introduced by von Mises? for the case of a buckled 
column. The procedure involves a set of successive approxima 
tions in dealing with the nonlinear terms of the governing equa 
tions; the first approximation has the form of the linear solution 
The illustrative problem presented here, a rectangular plate 


with all edges simply supported and under constant edge dis 
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stress. 


placements “uo and v% (see Fig. 1 inset), was chosen because pre- 
vious approximations (for example, references 3 and 4) and an 
exact solution® are available for comparison. The values of the 
effective width are in good agreement with those given by the 
exact solution (Fig. 1). 

For a plate with large deflections, the stress function F and the 
deflection w can be obtained from the von Karman equations: 


VIF = E(wry? — WrrWyy) (la) 
Viw = (h/D) (werFyy + WyyFrr — 2WryF2,) (1b) 


and the edge conditions for the problem stated above are 


w = (Oat all edges 
Wrz = Oatx = +(a/2) 
Wy = Oat y = +(b/2) 
Try += — Fy, = () at all edges 
e44/2 (2 
y = f _ [C1/E) (F, vF,,) —(1/2)w,?] dx = 
constant 
e45/2 
Uo = f [(1/E) (F,, vF,,,) (1/2)w,?] dy = 
b/2 } a 
constant 
in the usual notation. 
The linear solution, with indeterminate amplitude A, 
F = a(y?/2) (3a) 
w = A COS ux COS ny (3b) 


is the first approximation, where o is the average edge stress in 
the x-direction and 
B= mr/a / 
n = nx/b f 


Eq. (8b) is substituted in the nonlinear terms of Eq. (la), and the 
stress function F, obtained by solving the resulting linear equa 
tion, is 
F = o(y?/2) — (EA?/32) [(n?/u?) cos Qux 4 

(u?/n?) cos 2ny] (4) 


Expressions (3b) and (4) satisfy conditions (2). It may be noted 
that the constant % may be made equal to zero by the introduc- 


tion of an edge stress o, = F,, at the edges y = +(b/2); in the 
a/2 


present case the resultant oan dx = 0. The stress in the 
x-direction is (taken positive in tension): 

or = Fy, = o + (p?EA?2/8) cos 2ny (4a) 
and it may be seen that, as expected, the greater part of the load 
is taken by the portions of the plate close to the edges y = +(b/2) 
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To determine the amplitude A, Eqs. (4) and (3b) are subs 


tuted in the right-hand side of Eq. (1b) with the result 


View = —(hApt/D) }[(o/p?) + 
(EA?/16) (1 + 84) 


(EA/16) (84 cos 3ux cos ny + cos wx COs 3ny 


COS BX CO 


The solution, satisfying the boundary conditions of Eqs. (2 


terms of the order of A*is 


w = A cos ux cos ny + B cos 3ux cos ny 4 


where 
(A/h)? $/[3(1 v (1 + 67)2/(1 + £4) 
[(a/ der) 1] for ‘a 7 
B = —EhA*/(16D) [84/(9 + B?)?] 
C = —EhA/(16D) [1/(1 + 982)? ” 
— = [x*D/(hb?)] [8 + (1/8)]? 
D = Eh3/(12(1 — v*)];8 = w/n 


It is shown in reference 1 that for (a/b 1, the approximatio; 


38 = 1 can effectively be used; this value is therefore « mployed ir 
the remainder of this note 
The effective width }, is defined by the relation 


Or dy 


e414 
b.o(b/2) hh =h | 
J —b/2 


Use of Eqs. (4) and (6a) gives the first approximation as 

(b./b) = (1/2) [1 + (ocr/o4/2) g 
It is interesting to note, incidentally, that this expression is ident 
cal with the first two terms of a binomial expansion as obtait 
from von Karman’s effective-width formula*—i.e., 
+ [(oer/obs2) — 1] } 


2) [(ocr/or » — | tT » Sa 


(h, 5) = (o-,/or Jy -= +1 


1+ (1 


Eq. (8), however, is not sufficiently accurate except for stresses 


close to the critical (see Fig. 1), and the next approximation is 
therefore required. It was obtained by the procedure describ 

above, using in the right-hand side of Eq. (la) the expression for 
w given by Eqs. (6) rather than that of Eq. (3b). The new fur 


tion /, obtained after much algebraic manipulation, is then 


F = (oy?/2) + (EA?/25600) [di (cos 2ux + cos 2uy) 4 
dy cos 2ux cos 2uy + d3(cos 4 ux + cos 4uy) 4 


d,(cos 2ux cos 4uy + cos 2uy cos tux 


where 
d, = —800 — 12@ + S8la?/20) } 
dy = Ja - 9a? 25 
g 
d3; = 3a | 
| 
d,; = (48a — 9a?) /100 


in terms of the dimensionless amplitude parameter ? 


a=(1- v2)(A h)? ot 


The relationship between @ and o(4/2)ay, analogous to that giver 


by Eqs. (6a) and again obtained from Eq. (1b), is 


(oUF ar/@o =z i+ 0.375a(2 — 0.0225a@ tT 
0.091314a?2 — 0.00)9219a l 
where ov is the average stress along either unloaded edg 


fF 
I 


The final expression for the new approximation to the ¢ 


width is 


(b./b) = 1 av) (0.875a) (1 + 0.00059625a 10! 


— (Cer/ O(s 


The variation of the effective width with the average edge stress 
can now be obtained by solving simultaneously Eqs. (10a) and 
(10b); this was done numerically. The resulting values for tht 
effective width are plotted in Fig. 1 and are in good agreement 
with the exact solution given in reference 5 and with the appro™! 
mation given in reference 4 
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On Distributed Roughness as a Means of Fixing 
Transition at High Supersonic Speeds 


vf liam B F allis 


enior Aerodynamics Engineer, Convair, Ft. Worth, A Division o 
General Dynamics Corporation, Fort Worth, Texas 


ecember 13, 1954 


_ KNOWN WIND-TUNNEL TECHNIQUE for fixing the 


boundary-layer transition location at subsonic and low 
supersonic speeds has been developed through the use of distrib 
ited roughness. Recent Convair experiments in the Jet Pro 
pulsion Laboratory 20-in. supersonic wind tunnel at the Cali 
fornia Institute of Technology have shown that distributed 
roughness may be used effectively as a boundary-layer trip at 
igher Mach Numbers than generally anticipated.! 
rhe experiments were carried out on an approximately para 


lic body of revolution of fineness ratio 10. The nose of the 


model was contoured to a slight radius. Distributed roughness, 


pplied to the model nose, extended approximately 3!/2 per cent 
f the body length aft of the stagnation point. The roughness 
particles were aluminum oxide and were of grain sizes 320, 220 


150, 100, 80, and 60. A coating of “Trivar’’ baking varnish 
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Fic. 1. Variation of zero lift drag with roughness size 
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Fic. 2. Comparison of experimental and theoretical drag 


variation with Mach Number 


approximately 0.0005 in. in thickness, was used to fix the par 
ticles to the model surface 

Force data were obtained as a function of grain size at Mach 
Numbers 2.0, 2.5, 3.0, and 4.0 and for Reynolds Numbers in the 
range 6.80 X 10° to 10.5 & 10°; a minimum of three force meas 
urements were taken for each surface condition. Schlieren 
photographs were also obtained to enable a visual study of transi 
tion location 

The zero lift drag coefficient referenced to the frontal area of 
the body (3.142 sq. in.) is plotted against grain size in Fig. | 
As the roughness size was increased at a given Mach Number, 
the drag increased from a low value, corresponding to the ‘‘clean”’ 
model with partial laminar flow, to a constant value independent 
of grain size. The constant drag level indicates that the boundary 
layer was effectively fully turbulent from the nose aft and that 
the wave drag was insensible to changes in grain size. The grain 
size required to produce fully turbulent flow increased as the 
Mach Number increased. The upper boundary of the curves, 
where the drag again became a function of grain size, was not 
determined in the present tests 

The interpretation of the force data was further substantiated 
through examination of the schlieren photographs. An estimate 
was made of the body wave drag variation with Mach Number 
for the cases with fully turbulent boundary-layer flow, by taking 
the difference between the constant drag levels and the corre 
sponding theoretical turbulent skin-friction drags The open 
symbols in Fig. 1 represent the sum of the wave drag increment 
and the computed skin-friction drag, accounting for the location of 
transition as shown by the schlieren pictures. Satisfactory cor 
relation with the force measurements is evident 

A theoretical calculation of the body wave drag at Mach Num 
bers 1.5 to 3.0 was performed using Van Dyke's second-order 
theory and neglecting the slight nose radius. Fig. 2 shows the 
drag coeflicient plotted as a function of Mach Number rhe 
symbols represent force measurements obtained with transition 
induced at the roughness strip; the solid curve represents the ag 
gregate of the theoretical wave drag and skin-friction drag 
increments. It is concluded from the results that the distributed 
roughness as applied to the model provided a convenient and 
satisfactory method of producing fully turbulent boundary-layer 
flow throughout the Reynolds and Mach Number range of the 


tests. 
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Stagnation Temperature Control for a 
Blow-Down Supersonic Wind Tunnel 


T. A. d'Ews Thomson 
Department of Aeronautics, University of Sydney, Sydney, Australia 
December 14, 1954 


A METHOD OF MAINTAINING constant stagnation tempera 
ture during operation of an intermittent, blow-down, 
supersonic wind tunnel of short running time has been tested at 
the University of Sydney. Pillow' has treated theoretically 
the process that occurs when a fluid with nonuniform tempera- 
ture distribution passes through porous material. He assumed 
that the material was contained in a cylindrical, heat-insulated 
duct, that conditions were uniform over any of its cross sections, 
and that the conductive transfer of heat in both media was 
negligibly small. If the rate of heat transfer between the two 
media is sufficiently large that the temperature difference be 
tween them is negligible at any point, and the porous material is 
initially at uniform temperature, 7), then the temperature, 7>, of 
the emerging fluid remains constant and equal to 7), until 


m = M\(s/c) + (p/h)] (1) 


where m is the mass of fluid which has passed through the porous 
material, ¢ is the specific heat of the fluid, and p is its density; 
AJ is the mass of the porous material, s is its specific heat, and h is 
its apparent density. The temperature of further fluid leaving 
the material is 

Tm) = T\(m — M[(s/c) + (p/h)]) (2) 


if 7\(m) is the temperature of the entering fluid. 

The theory was used to design a unit for the 6- by 3-in. wind 
tunnel at the University of Sydney. The air was passed on its 
way from the high-pressure storage to the nozzle, through glass 
wool. This was packed between gauzes into a number of open- 
ended cells arranged, in honeycomb fashion, in a container. In 
let and outlet valves were operated to maintain constant both 
the rate of mass-flow through the container and the pressure in 
it. The temperatures of the air entering and leaving the glass 
wool were recorded and corrected approximately for some small 
fluctuations of pressure. 

The results of these tests (Fig. 1) show the effectiveness of the 
method. The outlet temperature was constant, within +0.5 
per cent, until a mass of air, 17|(s/c) + (p/h)], had passed through 
the glass wool, and this result was independent of the rate of mass- 
flow. 

With further flow, however, the results did not agree with the 
prediction of Eq. (2). The initial, sudden drop of temperature, 
due to the Joule-Thomson effect, which is conspicuous in the dis- 
tribution of inlet temperature, was not observed at the outlet. 
A number of factors were examined, which are neglected in refer- 
ence 1, and all were shown to have negligible influence on the 
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performance of the system tested except the effects of (a) no 


uniform packing density of the glass wool and (b) a possible erro, 


in the value of its specific heat. It appears likely that the qj 


crepancy between the experimental results and the prediction o 


Eq. (2) was due to a combination of these two effect 
Experiments and_ theoretical | 
even small amounts of leakage past the glass wool impair ser 
ously the performance of the system. Breaking up of the glas 
wool, causing choking and leakage, was noted after freque 
tests, and a more robust porous material appears desirable 
Initial uniformity of temperature of the porous 


important, and this can be obtained by the invers¢ 


material 1 


passing air of uniform temperature through the material betwe 


tunnel tests, 
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Aerodynamic Coefficients for Cambering and 
Aileron Motions 


E. F. Tyler 

Douglas Aircraft Company, Inc., El Segundo Division 
alif. 

November 30, 1954 


_ AERODYNAMIC FORCE COEFFICIENTS for an airfoil under 
going changes in camber, translation, and pitch (tw 
dimensional incompressible flow) were derived by Spielberg 
This note adds the extra coefficients that are necessary if th 
airfoil also has a control surface hinged about the nose. Th 
notation of Spielberg will be used. Let the airfoil be oscillating 
in the following motions: 


= ¢((1l — x? eto! 


which denotes parabolic change in camber of the mean chor 
with amplitude ¢ feet, measured at the mid-chord point, positiy 
down. 

S=U0 X<¢C 


= Bb(x —c ett x oC 


which denotes rotation of the control surface about its hinge | 
with angular amplitude 8 radians (relative to the main surface 
In the foregoing equations x is the dimensionless coordinate of 
point on the chord measured from the mid-chord and referred t 
b, the semichord. c = — cos ¢ is the value of x for the hinge lin 
The aerodynamic coefficient of the generalized cambering for 


due to control surface rotation is 


Np = —[C(k)/k?] Mi — (i/R)C(R)N2 — (1/k?) Ns — (1/R)Ms 4 
where 
N, = (1/7) (+ — ¢ + Sin ¢) 
No = (1/r) }[(a — ¢)/2] (1 + 2 cos ¢) 4 
(sin g/2 (2 COs ¢ 

N; = (2/37) sin? ¢ 
Ni = (1/r) [r — ¢ + (4/8) sin ¢ cos ¢ — (sin ¢ cos® ¢/9 
NV; = (1/7) }(3/4) (r — ¢) cos ¢ 4 

(sin ¢/2) [1 + (cos? ¢/2)| — (sim ¢/0 
The coefficient of control surface moment about its hinge line a 
to cambering is 
Ty = —[C(k)/k?]7T, — (i/k)C(R)T2 — (1/k?)T3 — (i/k 1 
where 
T, = (1/7) ((r — ¢) (1 — 2 cos ¢) + sin g(cos ¢ — 2 


estimates also showed tha; 


process of 
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The Lift of a Blowing Wing in a Parallel Stream* 


Heinrich B. Helmbold 
aff Engineer, University of Wichita, School of Engineering, 
Wichita, Kar 
December 13, 1954 
W \ JET IS ISSUED THROUGH a blowing slot in the down 
stream direction tangentially to the surface or at the trail 
edge of the wing as shown by Fig. 1, an additional lift, AL, 
result. Let p denote the air density, ) the span of wing and 
jet, h the thickness of the jet at the trailing edge, V the undis 


rbed velocity, and v the excess velocity in the jet 





Experience shows that the additional lift AL of a blowing wing 
is considerably higher, say 5 times, than the vertical component 
Al, = pbh{((V + v)? — V2] sin b: (1 
of the excess-momentum flux of the jet at the trailing edge (te 
This fact in itself says that an additional circulation must be pro- 
duced by the turning action of the basic flow on the jet. 
Integration of the pressure gradient in curved flow across the 
thin jet sheet in the manner of Fig. 2 gives the pressure increase 


p2 — fr = pl V + v+)*(h/r) (2 


f 
vhere V now denotes the local average velocity in the jet 
lhe pressure increase across a thin sheet of the basic flow with the 
same geometrical shape would be 


po’ — fi = pVr(h r) (3) 


] ( 
Thus an additional pressure increase across the thin sheet of the 


jet is produced by the excess of centrifugal acceleration in the jet 


po — po’ = p((V + v7)? — VJ (h/? (4 
lhe radial pressure distribution in the basic flow on the convex 
side of the jet is shifted by this amount with reference to the 
inalytical continuation of the radial pressure distribution in the 
basic flow on the concave side (see Fig. 2). Since the total 
pressure in the basic flow is constant in any case, the pressure in- 
crease of Eq. (4) is connected by Bernoulli's law to the velocities 


at opposite points of the thin sheet under consideration. 


V2?) /2] (5) 


bo — po’ 


se p2’ = pl[( V2" — 


his pressure increase may be interpreted as the pressure dis- 
ontinuity at a lifting surface where the velocity jumps from 
Veto Vo’. Then 

2 po — pro’ (V+? —-VWh 
a - — = (6) 


2 > hr + & Vo’ + Vo Zl 


is the lifting vorticity (circulation per unit length) of the discon- 
unuous surface that simulates the action of the thin jet on the 
basic flow. Because the static pressure is continuous through the 


1et 


t surface, the excess of total pressure in the jet is 


* + ; : 
his research was sponsored by the Air Branch of the Office of Naval 
Research under Contract No. N-onr 201(01) 


M 341 








Vv Ste 


Section of 





i blowing wing 








JET SHEET— 


Fic. 2. Jet sheet pressure and vorticity details 


Apt = pi ((V + 2)? — V?)]/2; i 
Hence, by Eq. (4), 
= 2Ap(h/r 
and by Eq. (6 

e = (4/p) [Api /( V2’ + V (h/r g 
The vorticity «€ must vanish at infinity downstream where the 
jet is completely turned in the direction of undisturbed velocity 
and where the radius of curvature has become infinite 


Now let x and y; denote the coordinates of a point on the vortex 
sheet substituted for the jet and assume that the slope dy;/dx and 


the curvature 1/r = d*y,;/dx* are small of the first order every- 
where. Then the vorticity € is small of the same order if Ap, 
(V2’ + V2) is of the zero order of magnitude In this order of 


approximation 
(Ve' + Vs)/2 = Ve 

if V., denotes the undisturbed velocity and 

e = 2(Ap:/pV. )h(d?y;/dx?* (10 
At the trailing edge the distribution ¢(x) is continuously connected 
with the distribution 7(x) of lifting vorticity on the airfoil, which 
means that the lifting vorticity y at the trailing edge does not 
vanish as with the airfoil with no jet. At the x-axis the whole 
vorticity distribution, y together with €, induces a stream slope 


w(x ] y(x’ e( x’ ; 
—_- — ——- . dx’ + . dx = 
J 2rt Jo* —% . x’ — x 
l € yx’) l 


Ap: d*y({x’') dx’ 
. dx + ea we (11 
2rI Jo* —# 2x pV.,?2/2 : dx’? x'—x 


where w,(x) denotes the local downwash velocity and c the chord 
of the airfoil. The boundary conditions are 


—w,(x)/V. = dy,(x)/dx (12 
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for 0 <x Scand 


= dy;(x)/dx (13) 


for x > ¢ In other words, at the airfoil, the flow must follow 
the slope of the airfoil mean line y,,(x), and, behind the airfoil, it 
must coincide with the slope of the jet mean line y,(x Elimina- 


tion of wj(x)/V., from Eqs. (11) to (13) yields two integral 
equations that determine the unknown functions +(x) and yj(x 
or 7(x) and ¢«(x), since «(x) is connected with y,(x) by Eq. (10). 

At the airfoil, the vorticity « induces an upwash angle Ae and 
an increase AV, of the x-component of the velocity, because with 
reference to the x-direction the jet is situated below the trailing 
edge. In order to compensate for the induced upwash angle, the 
original vorticity yo, with no jet, must be increased by an addi- 
tional vorticity Ay. 

It is supposed that the air quantity of an airfoil blowing jet is 
taken from the general flow and not from within the airfoil. 
Furthermore, it is assumed that a thin stream leaving it at the 
trailing edge will not influence the lift of the airfoil, provided the 
flow passing the interior of the airfoil experiences neither gain nor 
loss of energy; this is the case with no jet,v = 0. Witha blowing 
jet, the concern is only with the jet vertical excess-momentum 
flux AJ, at the trailing edge, which must be in equilibrium with 


the integral of the additional pressure forces of Eq. (4 
o= | pV. + AVr)e ds — (Al,/b) (14 


Here s is the curvilinear distance from the leading edge measured 
along the lifting surfaces, and AV, is the x-component of the in- 


duced velocity; hence 


Vo +t AV. = [(Ve2’ + V2)/2] (dx/ds 


The lift of the airfoil with blowing jet results from applying the 
momentum theorem to a control surface tightly enclosing the air- 
foil and crossing the jet normally at the trailing edge. 


L/b = [i a V., + AVe) (yo + Ay) ds +(AL,/b) (15 
From Eq. (14) this becomes 
L/b - ff, AV. + AVz) (yo + Sy) ds 4 
f pV. + AVr)eds (16) 


The lift of the airfoil with a blowing jet can also be determined by 
applying the momentum theorem to a control surface enclosing 
the wing and jet far from the wing everywhere and crossing the 
jet at a right angle (Kutta-Joukowski), 


L/b = pV [fp co tana + fi eas] (17) 
( P 
Since Eqs. (16) and (17) must be compatible, 
f, AV.(yo + Ay) ds 4 J AVzeds = 0 (18) 
0 


By subtracting the lift of the airfoil with no jet from the lift of 
the airfoil with blowing jet, we get the additional lift produced by 


AL/b = pV Lf; Ay-ds +4 / eds | 


and (18) this additional lift exceeds the 


the jet 


(19) 


According to Eqs. (14 
vertical component of the excess-momentum flux of the jet at the 


trailing edge by an amount 


AL — Al, 


} 
) 


=p | (V., + AVr)Ay + AVz-yo] ds (20 


0 


Therefore, the additional lift is greater than the downward com- 
ponent of the excess-momentum flux of the jet at the trailing edge, 
at least if the additional quantities AV, and Ay produced by the 
jet are positive everywhere along the chord. This will always 
be true with a positive (downward) deflection angle 6,, at the trail 
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ing edge, since then a positive vorticity ¢ results, which produ 

positive AV, and Ay along the whole chord. . 
With small deflection angles 6;¢ the induced velocity A] 

small of second order, and the additional lift is, to the first order 


of approximation, 


AL Al, 





Transonic Lift over the Front Surface of 4 
Two-Dimensional Diamond Airfoil 


Robert Wesley Truitt | 


Professor of Aeronautical Virginia F 
Institute, Blacksburg, Va. 


December 6, 1954 


_ THROAT STATION CONCEPT Of Tsien and Fejer,! combined 
with the Maccoll-Codd Method,? can be used to obtain th 


transonic chordwise lift distribution over the front surface of 


Engineering, 


symmetric diamond wedge at small angle of attack. TI 


Throat Station Concept is applicable to a diamond wedge of 


semiangle, 6, at small angle of attack, a, provided that, for 


the incompressible case, the fundamental zero-angle-of-attack 
flow character is unaltered. This condition is met approximately 
ifa/@ < 0.1 

The front surface pressure distributions of the diamond shoul 
follow closely the von Karman-Tsien rule up to and including 
Mach Number, J)... The 


is located just forward of the wedge shoulder 


the throat station critical throat 
station, (x/c)rs, 
The approximate front surface chordwise lift distribution is thus 
found easily in the Mach Number range (0 < JJ) < J/ 

By the Throat Station Concept, the local Mach Number dis 
will 


tribution, corresponding to the solution for /),,, remail 


virtually constant (frozen) in the subsonic range: 


{ , 2 
+ (Mo.,? 1)/[(y + 1)6]*/5} 


IA 


or 


c LQ. £h <. 


The Principle of Stationarity 
OM, 
Ol 
establishes the fact that the local Mach Number freeze persists 
into the supersonic regime 


Sor. = }(M,, — ] [(4 + 1)@)°'5} > 0 


where J/o,. ; is the upper limit of the freeze for the wedge sem 
angle, 6. 
Using Bryson’s approximation for the normal shock relatior 


1 — M,? Mi? — 1 


Moy. = M, is found fror 


The approximate rang 


the upper limit Mach Number, 


Eq. (4) by substituting 1., = 1, 
Mach Number freeze is thus 


é S bo S fou 


where |fo.r| = [fou | 

The front-surface pressure distributions are calculated throug 
the transonie range [Eq. (5)| by using the constant local Mat 
Number distribution in the Maccoll-Codd Method.* 





* This research was supported by the United States Air Force, thro 
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c. 


LIFT-CURVE SLOPE 


at 


Fic. 2. 


RESULTS 


rhe present method was applied to the theoretical incompres 
sible pressure distribution for the front surface of a symmetrical 
10) it | 


sure changes due to angle of attack 


diamond (@ = angle of attack rhe calculated pres 


ire shown plotted in Fig 
in transonic similarity form in the range 15 < & < 0.983 
rhe range of Mach Number freeze 


fy) < 0.988 Integration of Fig. | 


was found to bx ) O83 


gave the reduced lift-curve 
) 


slope shown plotted with other findings in Fig. 2 
CONCLUSIONS 


Based on the results of the present method, the following con 
clusions may be drawn 


1) The subsonic peak in the lift-curve slope is due primarily 


to the front wedge pressures. The subsonic reduced Mach 
Number at the peak (g = 0.983) corresponds to the attain 
ment of the Mach Number freeze 

2) As long as the Mach Number freeze persists, the curve of 
lift-curve slope versus Mach Number has a negative slope Phe 


present calculated negative slope at & = 0 is in excellent agree 


ment with the theoretical value 


l (Cr a d Cr a ail ‘= 2 y¥+1 


(3) The change with Mach Number of the pressure-drag co 


efficient due to lift 


Ca — Cay = (dCi/da) a? 


for a given angle of attack, has the same type of variation as the 
curve of lift-slope versus Mach Number 

(4) The second or supersonic peak in the lift-slope curve can, 
method, be anticipated for increasing Mach 
Wo, ,. Thus, both the 
sonic peaks are due primarily to the flow over the front wedge 


Willmarth’s 


» 


by the present 


Numbers above subsonic and super 


The abrupt drop in lift, at ¢ 1, shown in 
Fig , Was 
explained by him to be due to the supersonic flow over the rear 


experimental results for the complete diamond (se¢ 


half.4 This is seen to be in direct contrast with the present 
findings 
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On Determining Turbulent Boundary-Layer 
Separation in Incompressible and 
Compressible Flow* 


H. Schuh 

Division of Aeronautics 
Sweden 

November 92, 1954 


ckholm 


° 


Royal Institute of Technology, St 


I’ IS THE PURPOSE of the present paper to show that a new 
criterion of separation as found previously by the author for 


the incompressible turbulent boundary layer! can be applied also 


* The Dr. G. E. Gadd for communication of un 


published data and comments on experimental results 


author is indebted to 
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Fic. 1. Velocity ratio [m/u;(xo)| sep Causing separation of an 


incompressible turbulent boundary layer initially in a fully de- 
veloped turbulent state. 2, //(xo), and u,(x9) are the initial 
values of s [see Eq. (1)], the profile parameter //, and the velocity 
outside the boundary layer 1. (—k) is the mean gradient of 1. 


to the corresponding compressible case. This extension supports 
certain predictions of Lange,? Gadd,* and Gadd, Holder, and 
Regan,‘ but some new aspects are also obtained by the present 
considerations. First, the basic ideas of the incompressible case 
are briefly restated. A semiempirical method is used and 
separation is determined by the profile parameter // = 6)/é, 
where 6; is the displacement and 6 the momentum thickness 
The following differential equation is used for calculating //, 


written in the usual notation. 


(dH /dx) = A(H) BUE)Y (1) 


where 


Ms 


z= (u6./r) 


I = (2/u,) (du,/dx 


m is an exponent related to the skin friction; A and B are to be 
determined from experiments. So far the approach is similar to 
that in other methods;* * however, a new point of view is applied 
when using the available experimental material. Most of it 
was obtained on airfoils where laminar separation is followed by 
turbulent reattachment. Downstream of the separation point 
the profile parameter // falls quickly at first, later more slowly, 
before rising again under the influence of adverse pressure gra- 
dients."° In the region of rapid fall and some distance beyond, 
the turbulence of the boundary layer cannot be considered fully 
developed, and hence these regions are to be excluded in order 
to obtain a rational theory. This is supported by the following 
argument. Experiments on fully developed turbulent boundary 
layers at zero pressure gradients! do not fit with data from re- 
gions where the turbulence was suspected to be in a transitional 
state. Excluding such experimental data, the coefficient A in 
Eq. (1) becomes very small, and a solution could be obtained 
in a general form (see Fig. 1). zo, H(xo), and u(x») are initial 
values of z, H, and 1 in a fully developed turbulent boundary 
(—k) is the mean value of the gradient of u, which was 
The curve 


layer. 
used only in a term related to the small quantity A. 
for kzy = 
of the parameter kz justifies a posteriori the simplifying assump- 


corresponds to the case A = 0, and the small influence 
tion of a mean gradient in the term mentioned above. Separa- 
tion was assumed to occur at J] = ©, consistent with a vanishing 
shear friction at the wall. Agreement with available experi- 
mental results® !* !3 is good. According to Fig. 1, separation 
depends on the total increase of the velocity , from its initial 
value up to separation, but only in an extremely small degree on 
the gradient of 1; its influence disappears altogether for boundary 
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layers with a length of several feet and Reynolds Numbe; 
above 1.10%, when kz > 0.01. 
may possibly be limited to flow with monotonously increasing 


The new criterion for separatio, 


adverse pressures, since such were the conditions in the exper 
ments upon which the present method is based 

The effects of compressibility on boundary layers are of quant 
tative but not qualitative kind, if interaction with the flow oy 
side the boundary layer is disregarded. Hence it should } 
possible to predict separation of the compressible boundary layer 
by a diagram similar to Fig. 1. Donaldson and Lang: use 
Ap/q as a criterion for separation of the supersonic boundary 
layer, which is equivalent to using a critical velocity ratio as ; 
Fig. 1. Ap is the static pressure increase from undisturbed flo 
to separation or to a corresponding characteristic point, and 


is the initial free-stream dynamic pressure. Experiments , 


separation ahead of steps and concave corners (wedges) were us¢ 
as well as separation on plane walls due to incident shock waves 
since the separation point was not located in most experiments 
characteristic points in the static pressure distribution wer 
used as a criterion for separation. These points were the first 
peak in the experiments with steps and an inflection point in thos 
with wedges and shocks. When experimental values of A; 

were plotted against Mach Number, two curves were obtained 


for the two different criteria of separation.? In cases with inter 
action of the boundary layer with shock waves, Gadd, Holder 
and Regan? predicted that (Ap/qi)sEpP is only a function of Mach 
and Reynolds Numbers and is independent of the means by 


which separation is obtained 





In applying the criterion of reference 1 to the compressible 
case, it is noted that initial conditions in the experiments cor 
sidered were fully developed turbulent flow under constant pres- 
sure conditions, and hence the initial value of /7 depends only on 
Reynolds and Mach Numbers 
points of (Ap/qi)sEp against Mach Number can be expected to 


Consequently, the experimenta 


lie approximately on one curve for not too great a Reynolds Num 
ber range. In analogy to the incompressible case, the influenc 
of Reynolds Number on separation would be caused mainly, and, 
above a certain Reynolds Number (1.10"), entirely, by the Reyn 
olds Number variation of the profile parameter of the undis- 
turbed boundary layer, which is under constant pressure condi 
tions. The experiments of reference 4 indicate that the Reyn- | 
olds Number influence decreases with increasing Reynolds 
Number and vanishes at Re = 1.107. These predictions about 
separation of the compressible boundary layer are valid not only 


for shock-induced pressure rise but also for a continuous increase 
in pressure. 

Experiments with direct location of the separation point wert 
made by Fage and Sargent," Bogdonoff and Kepler, and Gadd 
Holder, and Regan.* The results showed that the critical pres 
sure ratios for incident shocks, wedges, and steps are independen 
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Fic. 2. Critical pressure coefficient (Ap/qi)sep for separatie 
of a compressible turbulent boundary layer initially in a fully 
developed turbulent state at constant pressure (undisturbed flow 
Valid for separation ahead of steps and concave corners and due 
to incident shock waves at plane walls. 
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f the m« yy which separation is obtained. From these ex- 


periments the pressure ratios referring to previously used char- 
in the pressure distribution can be correlated to 


icteristic points 
With the help of such a correlation, pres 


separation point. 
nts for separation were obtained from data quoted 


the 


sure coellcie 
‘1 references 2 and 17; they were plotted in Fig. 2 together with 
those from direct location of the separation point The point 


for incompressible flow in the figure was taken from Fig. 1 with 
H 1.35 (see reference 18) assuming initially constant pressure. 

The criterion of separation must be restricted at present to the 
three of separation to which the experiments of Fig. 2 re 


types 
fer. In one case of flow along a concave wall formed as a Prandtl- 
Mever expansion stream line, Drougge'® could detect no separa- 
tion at values of Ap/q = U = 1.81 


much higher than corresponding critical values in 


0.48 at and 2.76, and these 


values are 


t is not thought that this can be explained by the pres- 
A tentative ex- 


Fig. 2. I 

sure rise being continuous in these experiments 
lanation is that in this case a separated air layer cannot exist 
due to its interaction with the outside flow, and hence separation 

is prevented by extraneous forces This may possibly result in a 


type of turbulent flow of the boundary layer similar to those 
previously classified as not fully developed. Fig. 2 could also be 


ipplied approximately to bodies cf revolution,! but caution must 
be exercised for possible extraneous effects similar to those found 
Even so, Fig. 2 gives at least a lower limit 


on concave walls 


for the critical pressure rise 
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FORUM 345 


On the Effect of Reynolds Number upon the 
Peak Pressure-Rise Coefficient Associated 
with the Separation of a Turbulent Boundary 
Layer in Supersonic Flow 


E. S. Love 
Langley Aeronautical Laboratory, NACA, Langley Field, Va 
November 10, 1954 


ope PURPOSE OF THIS NOTE is to clarify the source of certain 
experimental information dealing with the peak pressure 
rise coefficient ?, associated with the separation of a turbulent 
boundary layer produced by forward facing st« ps and to present 


additional experimental data 





J ae eSern 

/0 /0 10 
Reynolds number 

Fic. 1. Variation of peak pressure-rise coefficient with Reynolds 


Number 


Reynold 
sf 1.55 


data at Vl = 55 that 


In a past investigation to determine the effect of 
Number on P,, the 


showed P, to be essentially invariant with Reynolds Number 


author obtained 


(See Fig. 1 Subsequently, similar results were observed at 


M = 3.03 and have been published in reference 1 In reference 


1 the data for \J = 1.55 were utilized in comparisons and were 


inadvertently attributed to the Langley 9-in. supersonic tunnel 
this has given rise to some question regarding the magnitude of 
the Reynolds Numbers involved (as high as 20 & 10°) in view of 
the maximum 
rhe data 


1.55 were obtained in a blowdown tunnel facility of th 


the limitation of 9-in. supersonic tunnel to a 


stagnation pressure of approximately 4 atmospheres. 
for Vf = 
Langley 9-in 


Supersonic Tunnel Section which is capable of 


stagnation pressures in excess of 30 atmospheres. Some addi 
tional data for higher Mach Numbers from tests at a later date 
in the same facility are presented in Fig. 1 These data also 
indicate only slight effect of Reynolds Number 
REFERENCI 
Lange, Roy H., Present Status of Information Rela to the P 
f Shock-Induced Boundary-Layer Sepa n, NACA TN 30¢ 14 


Transient Temperatures and Stresses in Plates 
Attained in High-Speed Flight 


). S. Przemieniecki 
Senior Technical Assistant, Structure »artment, The Br 
Aeroplane Co. Ltd., Filton, England 


November 9, 1954 


A SIMPLIFIED METHOD is presented for calculating transient 
temperatures and stresses in plates, such as cabin window 


panels or aircraft skins, 
All results are presented in the form of dimensionless 


subjected to a sudden aerodynamik 
heating. 
parameters so that these are applicable to different panel thick 


nesses and different materials 
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FIG Nomogram for local Stanton Number, .S? 


The solution of the differential equations of heat flow is ob- 
tained by the Laplace transformation method assuming zero 
heat transfer on the inside surface of the plate. The surface 
heat-transfer coefficient on the outside surface is determined 
from the aerodynamic conditions in the free stream just outside 
the boundary layer and the value of the local Stanton Number 
calculated from an empirical formula. 

The results for the sudden increase in kinetic temperature are 
generalized for the case of variable temperature by means of the 
Duhamel’s integral assuming constant surface heat-transfer 
coefficient during the heating process 


SYMBOLS 


v distance measured from the insulated surface, ft 
z distance measured from the middle surface, ft 
panel thickness, ft 


t time, sec 
1 T (x,t), temperature, °F 
oo. o kinetic temperature increase, measured from the initial tem 


perature, °F 


1 7 /¢o, dimensionless temperature parameter 

k thermal conductivity, B.t.u. per ft.? sec F./ft 

h surface heat-transfer coefficient, B.t.u. per ft.* sec. °F 

k thermal diffusivity, ft.* per sec 

] specific heat of air at constant pressure, B.t.u. per lb. °F 
KF local skin-friction coefficient 

p density of air, lbs. per cu_ft 

V free-stream velocity, ft. per sec 

n x 

\ « t/d?, Fourier modulus 

m k/hd, relative thermal surface resistance 

Py Prandtl Number 

I Reynolds Number 


Stanton Number 


| e—PNT, dN, Laplace transform 
J 0 


ith root of mw; cot w; 
on, ¢ thermal stress at ” and » 1, respectively, !t 
X a, (1 v) ‘aE, stress parameter 
x oi(l v) ‘akg, stress parameter at the heated surface 
I Young's modulus, Ibs. per sq.in 
a coefficient of thermal expansion, (in. per in.) per °] 


Poisson's ratio 


SOLUTION OF THE DIFFERENTIAL EQUATIONS OF Heat F 


The partial differential equations for one-dimensional hea 


flow in a panel facing a suddenly applied constant temperatyr 
gy expressed in a nondimensional form are:! 
0?7,/On? = OT, /ON 
oT, /on = 0 (atn = 0 9 
OT,,/On = (1 — T,,)/m (atn = 1 
T,(n, 0) = 0 
where T, = T/do, n = x/d, N = xt/d?, and m = k/hd. Th 


boundary conditions (2), (3), and (4) represent: 
(1) Zero heat transfer on the inside surface 


(2) Heat transfer on the outside surface controlled by the aero 


dynamic conditions in the free stream. This gives < 


condition ROT /Ox = h(@ — T), where the surface heat-transfer 


coefficient is given by 
h= Stepp V 
(3) Uniform temperature across the panel at ¢ = 0 (assume 


as zero for simplicity 
In Eq. (5) St is the local Stanton Number, c, is the specifi 


heat of air at temperature corresponding to the boundary-layer 


temperature near the heated surface, p is the density of air, and 


V is the free-stream v elocity just outside the boundary layer 


> 


The parameter m in Eq. (3) may be regarded as a measure of 


the relative surface resistance to temperature changes. For m = 


0 the temperature gradient at the heated surface is infinite, whict 
means that the outside surface reaches instantaneously its maxi 
mum temperature @, while for any other value of mm the temper 
ature gradient is finite 

The solution of the differential equations of heat flow can be 
obtained by the Laplace transformation method. Writing the 


Laplace transform as 


and applying the Laplace transformation to Eqs. (1), (2), and 


and solving the resulting ordinary differential equation for / 


have that 
ta = cosh (4 ‘pn p(cosh 1 Pp + mvp sinh yp 
To determine 7, from ¢, we use the Fourier-Mellin inversior 


theorem, which gives 


ae —~ Sil @; COS win \ 
7 =i1-—-—2 P . é - he) 
j Wy Tt SIN W; COS W, 
where w; are the roots of 
mw; = cot W; by 
If the relative thermal surface resistance mm = 0, then from Eq 
(9) cot w; = 0 or w; = qr/2 where g = 1,3, 5... and the expres 


sion for temperature distribution becomes 
- » 


| sin g7r/2 
i =ji— > COS G7rHn 2 r*g?N/4 


© q=1,3,5 q 


For other values of m roots of the transcendental equation mt 
cot w; can be determined from tables in reference 2 
The local Stanton Number (.S/) occurring in Eq. (5) depends 


1 ¢he 
1 the 


At 


upon the skin-friction coefficient, the Prandtl Number, an 
nature of the boundary layer, whether laminar or turbulent 
present, the existing knowledge on the Stanton Number 4l 
supersonic speeds is extremely limited, and the transition regio! 


is not easily defined at these speeds. It is, therefore, reasonad 


s a boundary 
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READERS’ 


(b) M=2-27 
2. Comparison of local Stanton Number given by Eq. (13 
with the experimental results of reference 6. 


assume that the boundary layer is completely turbulent, 


which gives a higher kinetic temperature rise, and use the sub 


sonic ¢ xpre ssion 


the von Karman relation 


St = (1/2)cer(Pr)— (11 


In the above expression cr is the local skin-friction coefficient, 
vhich is a function of the local Reynolds Number and the Mach 


the follow ing 


Number Experimental evidence suggests that 

formula is the most suitable for calculating cr (see reference 4 
0.472(1 — 1.12/logyw R)/(logy R)?*8 (1 + 0.2.17 r (i 

=|. ‘ 
































O2 


with nondimensional 


of stress parameter Y 


> Variation 


time parameter 


FORUM 34 


The Prandtl Number in Eq. (11) may be assumed as constant, 


since, according to existing data on the properties of air, 
from — di) to 


it only 


varies from 0.73 to 0.65 in the temperature range 


be important compared with possibk 


500°C., which may not 

errors from other sources In the correlation which follows, 
Pris assumed as 0.68. Hener 

St 0.268(1 — 1.12/logyw R logy, R S(1 + 0.2M 13 


local Reyn 


nomogram in 


The variation of the local Stanton Number with the 
olds Number and Mach number is plotted on a 
(13) gives good agreement with the recent experi 


Fig. | Eq 
results on boundary 


transfer in the turbulent 


These results® are plotted 


mental the heat 
layer on a flat plate in supersonic flow 


in Fig. 2 together with the corresponding curves calculated from 
Eq. (13). 


PHERMAL STRESSES 


It is shown by Timoshenko’ that thermal stresses in a flat plate 


at large distances from the edges, provided the plate is unre 


strained against rotation and expansion, can be calculated from 


Pe 
l4 

where 2 is measured from the middle surface of the plate. Sub 

stituting Eq. (8) into Eq. (14), noting that 7 = 7), and 


d {n — (1/2)], and integrating across the panel thickness, we have 


o,(l— vy . Sin & 
~ = ~ = —2 2, . 
ak@ ; @, > Sin Ww; COS & 
sin w I 
= COS n— 12(n - ° - 
cos & Sin &% ‘ m 
- P 15 
9 | 
where ,, is a dimensionless stress parameter. Putting » 1 in 


Eq. (15), we obtain an expression for the thermal stress on the 


heated surface 


o}( — sin &% 
Yo o I y. ~~ ; 
«\1 “ s 
ake ; a + SIN Wy; COS & 
$ sin w; 6 cos % 6 
T — COS & —_ f . 16 


aj 


The variation of this parameter with N is plotted in Fig. 3, while 
a typical stress distribution across the panel thickness for 1 /m 
38and N = 0.02and N = 0.2 is illustrated in Fig. 4 
VARIABLE KINETIC TEMPERATURE (1 
If the surface heat-transfer coefficient is assumed to be con 


stant during the heating process, the results for constant kineti 





COMPRESSION 2 
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TENSION 7 


A typical stress distribution across the panel thickness 
for N = 0.02 and N = 0.2 and 1/m = 3. 
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temperature increase ¢@p can also be extended for the case of vari 
able temperature ¢(t). Putting @) = 1 in Eq. (8) and using Du- 
hamel’s integral’ we have 

\ 
N,)] dM, 


ak. : . 
o( N,) _[F(n, N 
ON 


where 

: : ‘ W; COS WN 
Fin, N — N, - 
Sin W; COS W; 


(1S 


It should be noted that the variable temperature ¢ is expressed 
here in terms of the dimensionless time parameter. Furthermore, 
assuming that F(n, 0) = ¢(0) = O and integrating by parts, we 
have 

Op 


) dN (19 
ON, 


- N, 

which is more suitable for numerical computation than Eq. (17). 
Similarly, an expression relating thermal stresses and the 
time parameter V for the variable temperature can be obtained 


by substituting Eq. (19) into Eq. (14). 


ANALYSIS OF RESULTS 


Fig. 3 shows that the assumption of zero surface resistance 
(1/m = o) may lead to excessive conservatism in the panel de- 
sign. For example, 1/m = 3—an order of magnitude occurring 
in aerodynamic heating of glass panels—yields a maximum stress 
that is 19 per cent that of zero surface resistance. The point at 
which the maximum thermal stress occurs is a delayed function 
of time, the delay becoming increasingly pronounced as the surface 
resistance is increased. The maximum surface stress is an increas 
ing function of the parameter 1/m = hd/k. Thus, it increases 
with the surface heat-transfer coeflicient and panel thickness and 
decreases with the thermal conductivity of the panel material. 
This incidentally explains, for example, why thick glass panels 
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break more easily than thin panels under the same heating condj. 
tions 

The convergence of expressions for temperatures and stresses 
is slow if .V is small, so that in computing the stre parameter 
X’, a large number of terms had to be included for small values of 
N. The range of values plotted in Fig. 3 is especially suitable 
for low conductivity materials. For metallic materials the 
range is 1/m l and usually V > 0.2. (For a steel plate 0.25 in, 
thick V = 0.2 corresponds to 0.68 sec Thus it can be concluded 
from Fig. 3 that the thermal stresses in flat plates (unrestrained 
against rotation and expansion due to temperature changes) 
made of high conductivity materials are of relatiy ely small 
magnitude. 

For values of NV > 0.2 the analytical formulas for stresses and 
temperatures can be used directly, since retaining only the first 
two terms in the expansions results in an accuracy of about ] per 


cent, which is sufficient for most engineering calculations, 
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